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1 INTRODUCTION 

1.1 Problem Statement 

Spacecraft designers have been concerned since the 1960’s about the effects of meteoroid 
impacts on mission safety. Recent concerns have extended to the space debris environment, 
which typically displays more massive particles than the meteoroid environment for the same 
risk level. Additionally, the higher exposure area-time product of future space missions (e.g., 
Space Station) poses a more critical design problem than current short term missions. Finally, 
the inherent uncertainties in projectile mass, velocity, density, shape, and impact angle make 
the traditional deterministic design approach impractical. 

The engineering solution to this design problem has generally been to erect a bumper or 
shield placed outboard from the spacecraft wall to disrupt/deflect the incoming projectiles. This 
passive measure has resulted in significant structural weight savings relative to a single wall 
concept with the same protective capability. The problem, then, is how to efficiently design 
these protective structures so that the bumper disrupts the projectile without posing a lethality 
problem to the wall protecting the crew and equipment. 

Spacecraft designers have a number of tools at their disposal to aid in the design process. 
These include hypervelocity impact testing, analytic impact predictors, and .hydrodynamic 
codes. Perhaps the most widely accepted of these tools is impact testing, which has the advantage 
of providing actual spacecraft design verification. On the other hand, maximum test velocities 
are currently limited (8 km/sec) relative to maximum space debris (about 15 km/sec) and 
meteoroid (about 72 km/sec) velocities. Also, extensive testing is required to develop statis- 
tically significant trends for the large number of parameters associated with hypervelocity 
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impact Hydrodynamic code analysis can overcome the velocity limitation problem. However, 
this method is very computer (and time) intensive, and there is a fair amount of controversy 
involved in the selection of appropriate codes and code-specific par amet ers. 

Analytic impact predictors generally provide the best quick-look estimate of design 
tradeoffs. Their use is constrained by the limitations of the testing from which they are 
experimentally derived, the assumptions used in their theoretical derivation, or the regression 
analysis used in their statistical formation. However, analytic predictors may provide infor- 
mation that is clearer than that obtained from the examination of experimental results. 

The most complete way to determine the characteristics of an analytic impact predictor 
is through (nonlinear) optimization of the protective structures design probl em f o rmu lated with 
the predictor of interest. Optimization techniques provide analytic or nume ric al solutions 
depending on the nature of the predictor, the problem formulation, and the technique used. 

1.2 Contract Pu rpo se 

The purpose of this contract is to provide Space Station FREEDOM protective structures 
design insight through the coupling of design/material requirements, hypervelocity impact 
phenomenology, meteoroid and space debris environment sensitivities, optimization techniques 
and operations research strategies, and mission scenarios. Major findings from contract 
inception to the beginning of this study are detailed in References 100-105 and are shown below: 
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Posynomial Regression Can Be Performed To a Statistically Significant Level for 
Hypetvelodty Impact Test Databases, i , . : v > 
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13 Study Goals 

The goals of this study are to: 



e to the *new" space 


L Perform Space Station protective structures design 

debris environment definition, ;■ f" yv-.' ", . , , K A.. /: -.' ', 

2. Incotporate the unique methodology developed into a user-friendly, menu-driven PC tool. 




for Space Station protective structures designers, 

4. Assess the hypervelocity impact test samples from a daraage/penetration standpoint 

5. Analyze projectile shape effects on protective structures design. > , • 


The period of performance for this effort is 2-28-90 through 6 -30-91 


e Scope 


6. Perform preliminary advanced shielding development work in the area of multiple bumper 

configurations., v - " , t' \ ^ ,, y- , ' 

7. Develop discrete protective structures design optimization methods. 


1.4 Study Approach 

The methodology presented in this study is sufficiently general for application to various 
spacecraft configurations and impact environments. The baseline scenario investigated is for 
the Space Station Core Module Configuration and space debris environment with the following 
specifications: 5 % space debris growth rate; Space Station operation period fro m 1995-2004; 
460 km Space Station altitude; 28.5 degree Space Station inclination; 0.97 total Core Module 
Configuration probability of no penetration; 588 m 2 total Core Module Configuration debris 
area; 10 cm bumper/wall separation; 0 degree impact angle (normal); 6061-T6 aluminum alloy 
bumper, 2219-T87 aluminum alloy wall; and 9 km/sec average impact velocity. 

Because other approaches involve the analysis of existing protective structures designs, 
the design methodology presented here is unique. The process begins with the definition of the 
space debris environment to determine the critical design projectile diameter and density. The 
design problem is then formulated in terms of a hypervelocity impact predictor as a weight 
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minimization function of the independent (or designer controllable) variables. These variables 
gener all y include bumper/wall material properties and thicknesses. The protective structures 
system is then globally optimized using the Geometric Programming technique. Sensitivity 
analyses are performed to investigate the effect of changes in the system parameters on the 
optimal design. Several hypervelocity impact predictors are analyzed, including the Wilkinson, 
Burch, PEN4, and Nysmith models, as well as combinations of these models. 

1.5 Study Results 


1. Goal 1 was completed using technology-specific optimization techniques. Results are 

given in Section 2. " ! - , . 

2. The user-friendly, menu-driven PC tool of Goal % is called PSDOC (Protective Structures 

Dc$ign Optimization Code) and was delivered with 4o<riu^n^tiOii in August 1990. Several 
updated versions have been delivered in the interim. An overview of this tool is presented 
in Section 3. ' ' ' . ' 

3. The Monte Carlo simulation tool of Goal 3 has: been planned and is currently under 

development. This i$ discussed in Section 4. ■ . . ^ . 

4. Hypervelocity impact, samples (Goal 4) have been evaluated as discussed in Section 7. 
$AIChas also performed additional posyhomial regression analyses on hypervelocity impact 
test data which has been delivered in a White Paper andas part of a Dissertation. gig ? 

5. SAIC has developed appropriate regression and optimization tools to satisfy Goal 5. This 

is presented in Section 8. y* / .. '' ' ' 4, 

6. The development of advanced shielding concepts ha$ been performed (Goal 6), and will 

continue to be performed as part of future work on this contract. Section 5 includes results 
from this effort x , ' + # 

7. The development and applications of discrete protective structures design optimization 

techniques (Goal 7) is complete and is presented in Section 6. " : 
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1.6 Major Findings of This Study 


1. Globa! analytic nonlinear design optimization can be performed for the projectile melt/Vaporization region 
(Wilkinson), for normal impacts in the projectile shatter region (Burch), and for the Nysmith predictor using 
Geometric Programming. '*• ' ' ' ' 

2. For the predictors investigated, the optima! ratio of bumper mass per unit area to total mass per unit area 
may vary with mission, environment, projectile mass, and velocity regime. • ~ . - ' " ' 

3. There is a large incwByi’'ifor increasing the bumper/wall separation from 10 to 15 cm for ail predictors 
investigated - ' ... . > , ' r. = 

4 , Ail predictors reflect lwimasing design sensitivity to projectile #umeter and decreasing design sensitivity 
to bumper/wall separation. J U ... 

S. The Wilkinson and Nysmith predictors reflect increasing design Sensitivity to projectile velocity, while the 
Burch relationship is decreasing. 

6, For the combined predictors, 2Q11-T8 is the preferable aluminum alloy bumper choice for the baseline 
parameters, .... | H w - ' ' , ' ' 

7, For the combined predictors, increasing the bumper/wall separation from 10 to 15 cm reduces the minimum 
module weight by 25%. (. . . | ' ' ' - 

8 . Minimum CMC weight is very sensitive to space debris growth rate above 7% and Space Station altitude 
below 1000 km for the combined predictors. ' ... .7 ; I 

9,. .CMC pfoiictiyl^ greatly on mission duration for the combined predictors. 

10. For the combined predictors, increasing the CMC mission risk from 3% to 5% reduces the minimum module 
weight by about 30%. ' j - - ... Ilifi 

rill: Global (and sometimes analytic) optimizatibh of discrete posynomial programs can be performed using 
dual approaches coupled with partial invariance techniques. .. ' ' . e 

12. Primal methods require less "pencil and paper" effort than dual methods and are more easily applied to 
most problems, 1118 g llljitll | ' , . ,.Z...7 

13. Primal methods do not generally obtain global solutions for the discrete posynomial program. 

14, The dual method may bb advantageous in cases where the objective function may be sufficiently separable, 
since posy separable programs do not require solutions of coupled nonlinear equations in the dual-to-primal 
variable transformation. ,,,, , , ~ ui “f - r J 

15. The Monte Carlo simulation tool is feasible from a development standpoint and appears to have advantages 
oyer current expected value models. II . •• 

16. Aiii approximation to a nonstationary Poisson arrival process for impact events spears to be sufficient 

17, Both the Wilkinson and ballistic PEN4 predictors may be extended to multiple bumper models. 1111 

1$, The multiple bumper Wititihson predictor optimization problem is a 0 degree of difficulty posynomial 
programming formulation. * •• . • 
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2 ANALYSIS OF NEW SPACE DEBRIS ENVIRONMENTS 
2.1 Earth Orbital Space Debris and Meteoroid Environs 

The space debris environment model chosen for this study is due to Kessler 75 . The major 
dependencies considered involve space debris growth rate, spacecraft operational period, 
mission altitude and inclination, spacecraft debris area, orientation, and probability of no 
penetration. 

The space debris flux is given by Kessler as 

F(D,h,i,t,s) = Bty(h,s)y(i) (F,(D )g l (t) + F 2 (D )g 2 {t )) [ 1 ] 


where 


<KM) = sy(<l>i(M)+l) 

i /■ » v in (A/200-*/140-1.5) 

^(h,s)= 10 ' 


F,(D)=1.05(10' S )/Z) z5 


FJD ) = 7.0(10 lo )/(D + 700) 6 


£ i (0 = (1 +2 P )' -1985 

g 2 (r) = (l+F )'- 1985 

The spacecraft inclination factor for 28.5 degrees is 0.9135. 
The cumulative flux N is given by 

N= [ T FAdt 

Jo 

which may be approximated using one year intervals by 


N=A 't F(D,h,i,t,s(t)) 



[2] 

[3] 

[4] 

[5] 

[ 6 ] 

[7] 

[8] 
[9] 
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A Poisson arrival rate for space debris gives 

/><, = «■" [ 10 ] 

A closed form solution for D may be accurately found for particle diameters much smaller than 
700 cm. This is given by 


D = 


1.05(10" 5 )(G,) 


W 


-5.9499(10- 7 )G 2 - 


ABWtj 


where 


for j= 1,2. 


The average projectile mass density is given in gm/cm 3 by Kessler as 

p p = 2.8 for D £ 1cm 


p = 2.8/D 074 for D > 1cm 


[ 11 ] 


[ 12 ] 

[13] 

[14] 


This relationship is shown in Figure 2.1-1. 


For an orbital inclination of 28.5 degrees, the non-normalized impact velocity distribution 
is given by 

/(V) = (14.46V- V 2 )(i8^^' 18W1 ^ + 0.67e^ _9 - M5y3 925)l )+0.01 16(28.91 V - V 2 ) [15] 
The normalized impact velocity distribution is given by 

[16] 

Jo f(YW 
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This distribution is shown in Figure 2.1-2 for i = 28.5 degrees. Finally, the impact angle is 
given as a function of impact velocity as 

0 = cos” l (— V715.4) [17] 

This relationship is shown (with uncertainty bounds) in Figure 2.1-3 for a surface parallel to 
the CMC velocity vector. 



Figure 2.1-1. Space Debris Particle Density vs Diameter 
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Figure 2.1-3. Projectile Impact Angle From Normal of Surface Oriented Parallel to 

CMC Velocity Vector vs Impact Velocity 
The total meteoroid environment flux-mass model is given by Cour-Palais 32 as 

Log 10 (JV,) = -14.339 - 1 .584 Log 10 (m ) - 0.063(Log 10 (m )) 2 [18] 


for 


m € [ 10 ~ 12 , 10 " 6 ] 


and 


L°gio(ty) = -14.37 - 1.213Log 10 (m) 


[19] 
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with shielding factor 


m e [10^, 1] 


„ l+cos(<|>) 

n — 


sin(<)>) = 


R 


R+h 


[ 20 ] 

[ 21 ] 


and gravitational defocussing factor 

043 

G=-— + 0.57 [22] 

N r 

where R is the radius of the shielding body (= 6378 km for Earth), and N R is the spacecraft range 
from the Earth’s center in Earth radii. The velocity probability distribution for meteoroids is 
shown in Figure 2.1-4. The average mass density is 0.5 gm/cm 5 , with average particle velocity 
of 20 km/sec. 



Figure 2.1-4. Meteoroid Velocity Probability Distribution 
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2.2 Measures of Design Effectiveness 

The traditional measure of protective structures design effectiveness is the probability of 
no penetration of the pressure wall. This measure generally accounts for the risk associated 
with the particle size, impact velocity, and impact angle. It may also include spall factors to 
account for impact scenarios where penetration does not occur, but spallation does. 

The probability of no penetration of a protective structure generally does not reflect uncertainties 
in the environment, in particular, the particle shape, density, and diameter. These uncertainties 
may be estimated by establishing confidence intervals about the expected probability of no 
penetration. 

2.3 Potential Protective Structures Design Approaches 

Active Design 

Active design includes debris mitigation and removal. Debris mitigation is the design of 
spacecraft and launch vehicles to minimize the amount of debris generated through operations. 
Debris removal includes the entrapment and/or possible destruction or disposal of debris. 
Passive Design 

Passive protective structures design is the placement of shields permanently spaced 
outboard from the pressure wall to disrupt the incoming particle. One approach to providing 
design insight is through the use of Geometric Programming (GP). 

GP is a particular nonlinear programming (NLP) technique formalized by Duffin, Pet- 
erson, and Zener 42 in 1967. It is practiced by engineers, scientists, and mathematicians alike. 
To appreciate the elements of GP requires a short mathematical presentation. 

The prototype Geometric Programming problem is formulated in terms of posynomials - 
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polynomials with positive coefficients, positive-valued independent variables, and real expo- 
nents. The problem is to . 


min/= X Ci IT x? 

i-J y- 1 


[23] 


subject to 



Obviously, this is a great restriction in applicability, since not all NLP problems may be for- 
mulated in terms of [23] and [24]. For problems of this form, including nonconvex programming 
problems, GP provides the globally optimal solution. 


One approach to solving this problem is to consider the dual problem, as justified by the 
Arithmetic -Geometric Inequality. The dual Geometric Programming problem is given by 


maxv(5)= FI 

i«l 


Y< 


(finrme?)'')) 

'=» 5;, 


[25] 


with 



IS, = 1 [27] 

i* 1 ‘ ■ 


ft- I*, / = l,2,...,p [28] 

j-1 ■ -r - - 

Qeariy, equations [26]-[28] represent k+p+1 equations in n+p+m 1 +m 2 +...+m p unknowns. If 

k + l> n+ & m, [29] 
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then the system is overspecified. If, in addition, the system is inconsistent, then the problem 
formulation or model selection must be reconsidered. If 


k + l= n+ £ m, 


[30] 


/-i 


and the system has nontrivial determinant, then a unique solution for the dual variables exists. 
If 


k + l<n+ £ m , 


[31] 


t-i 


then the system is underspecified. The Geometric Programming degree of difficulty is given 
by 


DOD =n—k—l+ £ m, 


i - i 


[32] 


Optimal dual variables for systems with positive degree of difficulty may be found by using a 
number of techniques, including search methods. Once the optimal dual variables are deter- 
mined, they must be converted back to the primal variables using the relationships 

/o = v(5*) [33] 


c, IT x? = 8‘/o i = 1,2, ...,n 


j - 1 


o'* 


li,c a Ilx/ = 8 a / = l,2,...,p 


[34] 


[35] 


Note that this dual-to-primal conversion involves n+p nonlinear equations, and therefore rep- 
resents a potentially difficult problem to solve in its own right 

Now, if the number of terms in the objective function (n) is large, and the number of 
independent variables (k) is small, a large degree of difficulty problem often ensues (particularly 
in a problem with few constraints). In these cases, solution of the dual problem may be quite 
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lengthy, and a primal method may be in order. This strategy is further justified when gradient 
methods are used, because the first and second (and higher-ordered) partial derivatives of the 


independent variables are easily pven as: 


dxi 



[36] 


^Xca^cv,)-' ru? 


[37] 


Based on the relatively large number of recent applied Geometric Programming articles, 
it is apparent that GP possesses a fairly high utility, particularly in the area of structural design. 
Because GP is the only NLP technique which offers the guarantee of a globally optimal solution 
for certain nonconvex problems, it should be considered more widely in practice. Additionally, 
for zero degree of difficulty problems, GP can provide an analytic optimal solution for the 
objective function and independent variables. This attribute provides greater insight for the 
system designer than that obtainable by other NLP techniques. Finally, the values of the dual 
variables may provide very crucial design information alone in terms of the physical parameters 
of the problem at hand. 

Since its inception, GP has been widely applied to structural design optimization problems. 
These problems may involve dynamic and static loadings, both determinate and indeterminate. 

"i’sz::. . 

The posynomial property of weight minimization for structural design problems matches nicely 
with the GP technique. Additionally, since many structural design optimization problems 
include a large number of independent variables, this reduces the degree of difficulty for the 
GP process (see equation [32]). 
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Recently, GP has been found to be widely applicable to the optimization of spacecraft 
protective structures using analytic hypervelocity impact models. 100-105 The posynomial nature 
of these predictors is not unusual, since many physical phenomena may be attributed to a 
geometric model. 

The basic optimization problem is a weight minimization problem of the protective 
structures. It has been shown 105 that for spacecraft structures with low curvature and relatively 
large diameter, it is sufficient to minimize the total mass per unit area given by 

W=i p,t, [38] 

i«l 


In particular, this is true for the Space Station Core Module Configuration. Increasing the 
complexity of the weight objective function by accounting for specific configurations only 
serves to increase the complexity of the optimization technique and convergence time unnec- 
essarily. No improvement in accuracy is achieved. 

Three hypervelocity impact predictors, developed in the 1960’s and displaying different 
attributes of Geometric Programming are due to Wilkinson 160 , Burch 29 and Nysmith. 115 

The Wilkinson predictor is a piecewise differentiable model given by 


0.364P 3 p,ycos(9) 
M 2 P2 


for ^£1, 
Pl'l 


[39] 


O.364£) 4 p 2 Vcos(0) „ 
t- —S. for 

Pl^lP2 P 1 * 1 


DP '>1. 


[40] 


Under condition [40], the dual Geometric Programming objective function is given by 


v(8) = (p 1 /5 1 ) S, (c,/5 2 )‘ l 


[41] 


0.364D 4 p 2 y cos(6) 
M 2 Pi 



[ 42 ] 
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5j + §2 = 1 

[43] 

5,-S J = 0 

[44] 

Equations [43] and [44] together imply 


8, = = 1/2 

[45] 

The minimum weight and globally optimal thicknesses are given by 


1.207D‘p,(^f 
W » m 

[46] 

0.604 D’pJfrY 

Sp, 

[47] 

0.604D 2 p,( 1 * * * V '^ (8> ) 

* Sfc- 

[48] 


Thus, the globally optimal algorithm for the Wilkinson Predictor is 


1 Determine r t and t^from equations {47] and [48]. 


2, Compute 


3. If |f|p > I, then quit The optimal design is (f, , t). 


S 1, the optimal design is Hv^pTF 


Figures 2.3-1, 2, and 3 show the optimal design values of minimum system mass per 
unit area, and optimal bumper and wall thicknesses vs projectile diameter, bumper/wall sepa- 
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ration, and projectile velocity, respectively, for the Wilkinson predictor. In Figure 2.3-1, the 
projectile density varies with diameter according to equations [13] and [14]. In Figure 23-3, 
the impact angle remains constant at 0 degrees (normal). The optimal bumper and wall 
thicknesses for the Wilkinson predictor are approximately equal due to the similarity in bumper 
and wall material densities (see equations [47] and [48]). 



Projectile Diameter (cm) 


Optimal Bumper Optimal Wall Minimum System Mass 

Thickness (cm) Thickness (cm) Per Unit Area (gm/cm2) 


Figure 23-1. Optimal Design Value vs Projectile Diameter for Wilkinson Predictor 
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o 2 4 6 8 10 12 14 16 

Projectile Velocity (km/sec) 

Optimal Bumper Optimal Wall Minimum System Mass 

Thickness (cm) Thickness (cm) Per Unit Area (gm/cm2) 


Figure 2.3-3. Optimal Design Value vs Projectile Velocity for Wilkinson Predictor 


The normal impact predictor for the Burch model is given in functional form as 

(snsr 

h s 071 

where 

F, = 2.42(r,/£> )"* 33 + 4.26 (r,/D ) 033 - 4. 1 8 
Equation [50] may be approximated by 
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K = Fi 71 = 2.8 (/,/D l 0 - 57 + 1.58(r,/D)^ 
Then W is given in posynomial form as 

W = p L /j + p £K 


where 



The dual Geometric Programming problem is to maximize 


v(5) = (pA) 


r 2 . 8 p 

r 1.58p 2 CD 0J7> 

S* J 

83 J 


subject to 


6 t + 0.575J -0.5763 = 0 


S 5, = 1 

i-i 


Equations [55] and [56] may be partially solved to give 

52 = 2.33(1-1.5753) 

5, = 1.33(253-1) 

Since the dual variables must all be positive, we have 

0.5 < 63 < 0.64 

Thus, the one degree of difficulty algorithm is given by: 


[51] 

[52] 

[53] 

[54] 

[55] 

[56] 

[57] 

[58] 

[59] 
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1, Vary 83 from 0.5 to 0.64 to find the max v(8). 

2, Using the corresponding 5,, solve for 8, and 8*. ' 




3- ^i # -5jV 0 (8yp,. 

J| | | j|| v 0 (8)-p 1 f 1 , . 

U . ; ■ ^ '•“-'pa 


Figures 2.3-4, 5, and 6 show the optimal design values of minimum system mass per 
unit area, and optimal bumper and wall thicknesses vs projectile diameter, bumper/wall sepa- 
ration, and projectile velocity, respectively, for the Burch predictor. Figure 23-4 reflects a 
constant projectile density as given in equation [13]. In Figure 23-6, the impact angle remains 
constant at 0 degrees (normal). 
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Projectile Diameter (cm) > 

Optimal Bumper Optimal Wall Minimum System Mass 

Thickness (cm) Thickness (cm) Per Unit Area (gm/cm2) 


Figure 2.3-4. Optimal Design Value vs Projectile Diameter for Burch Predictor 








Bumper/Wall Separation (cm) 

Optimal Bumper Optimal Wall Minimum System Mass 

Thickness (cm) Thickness (cm) Per Unit Area (gm/cm2) 


Figure 2.3-5. Optimal Design Value vs Bumper/Wall Separation for Burch Predictor 
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Figure 2.3-6. Optimal Design Value vs Projectile Velocity for Burch Predictor 
The Nysmith equation was developed for meteoroid impacts and may be written 


with inequality constraints 
and 


5.08V OJ78 D z ” 

^0.528^1.39 


£ s0 - 5 


[60] 


[61] 
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Substituting equation [60] into [38] results in 

5.08 V 0278 # 2 - 92 

w = t + 

,f 1 jOJ2*^1.39 

The problem constraints may be rewritten 


[62] 


[63] 


[64] 


^ 21.72V a527 £> 3636 

*1 ^ £ 2.633 


[65] 


The first step in this analysis is to determine when the problem is feasible. This corresponds 
to the question: When is the constraint set defined by [64] and [65] nonempty? Clearly, this 
is the case if 


D ^ 2 1. 72 V * 327 /} 3 636 
2 * 5 2-633 


[ 66 ] 


or 


A more usable form is given by 


D <, 


0.2395 

y °-2 


5 2>4.184DV 02 


[67] 


[ 68 ] 


The conditions of existence of a local (and thus global) optimal solution to the problem 
will now be established. 

If 


D ^ 0.235 V^ 2 


[69] 
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then the optimal solution to the problem exists and is given by 


1.907V 0182 D 19 * 

[70] 

*1* — j0.9I 

3.613V 0 ln D 191 

[71] 

% — j0.91 

5.520V OI82 D 1 - 91 

[72] 

Wi- s «. 


Note that the ratio of optimal bumper thickness to total thickness is 0.345. The corre- 
sponding ratio for the wall is 0.655. Thus, provided the values of the systemic parameters satisfy 
[69], these ratios are constant." 

Finally, notice that we provide optimality conditions for most of the feasibility region. In 
fact, it is now only necessary to determine the existence of optimal solutions in the interval 

0.235 V~* 2 <,D<, 0.245 V -02 [73] 

Figures 2J-7, 8, and 9 show the optimal design values of minimum system mass per 
unit area, and optimal bumper and wall thicknesses vs projectile diameter, bumper/wall sepa- 
ration, and projectile velocity, respectively, for the Nysmith predictor. Figure 23-7 reflects a 
constant meteoroid density. In Figure 2.3-9, the impact angle remains constant at 0 degrees 
(normal). 
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Projectile Diameter (cm) 


Optimal Bumper Optimal Wall Minimum Total System 

Thickness (cm) Thickness (cm) Thickness (cm) 


Figure 2.3-7. Optimal Design Value vs Projectile Diameter for Nysmith Predictor 
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Figure 2.3-9. Optimal Design Value vs Projectile Velocity for Nysmith Predictor 


We now consider the combination of impact predictors corresponding to ballistic, pro- 
jectile shatter, and projectile melt/vaporization regions. The optimization problem is first 
formulated and then solved for these three impact regions. These optimal solutions are then 
integrated into an overall optimal solution. The predictor equations chosen are based on previous 
work performed by Boeing. The ballistic, projectile shatter, and projectile melt/vaporization 
predictors are given by the PEN4, Burch, and Wilkinson models respectively. 
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The PEN4 model in functional form is given by the following set of equations: 


c,p p Y 31 f 0.28 ID p Y* /m 

{ * J 

[74] 

a —b 

[75] 

Cl= 7+d 

a = 1.33V a *J>J 

[76] 

b = 85, i I,£- 3 '” <,o> /cos(0) 

[77] 

c = l33R 2 pPp 

[78] 

d = R p t lPl / cos(0) 

[79] 

This set of equations is valid for 


V <V / +4000 

[80] 

where 


V 7 =4100 if r,/D <0.4 

[81] 

V f = 4986(r,/D ) 021 if r,/D >0.4 

[82] 

When equations [74]-[79] are substituted into equation [38], a one- 

dimensional search is per- 

formed on t! with initial point 


3.i2jpo“ty- 

h = 0.16625V 2 /?;p, cos(e)— 

[83] 


corresponding to t 2 = 0. When a local optimal solution is determined, condition [80] is checked 
to determine if the ballistic region is appropriate for consideration. 
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The Burch model is actually two separate predictors, one for normal impacts, and one for 
oblique impacts. The normal impact predictor is given in functional form in equations [49] 
through [53]. The oblique Burch predictor is formulated in terms of flight path and normal path 
penetration as 


( F,+ 0.63^2 'j (C/V) 16T7 (D/5) V7 

l N f J 

[84] 

where F, is as defined in [50] and 


F 2 = 0.5 - 1.87^/D ) + (5r,/Z> - 1.6)% 3 + (1.7 - 12r,/D )x 

[85] 

X = tan(0) - 0.5 

[86] 

The weight minimization problem may then be formulated as 


W = PA + M 

[87] 

subject to 


N n <, 0.85 

[88] 

where 


N, = F 3 (D/I 2 )(C/Vf 

[89] 

F 3 = 0.32(r,/D f 6 + 0.48(r,/D) ,/3 sin 3 (0) 

[90] 


and tj is given by [84]. This problem is solved using an exterior penalty function technique 
with objective function 

$(0 = W + 5 K(N n - 0.85) 2 [91] 

where 

8=1 if N n -0.S5>0 [92] 
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8 = 0 if N n - 0.85 <0 [93] 

A random search with a 99% confidence interval of 0.01 inches is performed, and K is increased 
until 

8K(N n -0.&5) 2 <>E [94] 

The random search interval for tj is specified by using the single plate equation 

t l =K l m°- 35 Y 6 V ojn5 [95] 

0.816 


K — ' 

1 ^ 1 / 18 ^ 1/2 
e p. 


[96] 


The interval is then given by [0,t,]. 

Due to the discontinuities existing between the three impact predictors, an integrating 
algorithm must be developed. This algorithm is included for fixed velocities. 


1 . Compute optimal design for PEN4 predictor, 

- : 2. Check against PEN4 constraint [80]. 

v 3, If satisfied, die optimal design is 

4. Otherwise, compute optimal designs for Burch 
o _ and Wilkinson pi|(iictors, and (*iy * 2 y) respectively. 

^5. Compute Wilkinson wall induced by optimal Burch bumper, 

- 6. Compute Burch wall induced by optimal Wilkinson bumper, 

|3g|§3 
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Once the optimal bumper and wall thicknesses are determined for each velocity, the 

integrated optimal bumper and wall thicknesses are found from 

. v 




tJVMVMtVW for i= 1,2. 


[97] 


Real Time/Reactive Design 

Real time and reactive protective structures design refers to the concept of performing 
design in orbit through the use of smart structures, smart materials, or the combination of passive 
and active design techniques. The real time design approaches may be accomplished through 
particle sensing either before or during impact. Impact particle mass, velocity, angle, and 
location prediction is performed to provide the necessary algorithmic information to the 
structure/material controller. The material/structure is then configured to defeat the specific 
impact scenario anticipated. Real time/reactive protective structures design provides the most 
flexible and safest design alternative available, but also stresses technology the most. 

2.4 Aluminum Alloy Bumper Materials 

A comparison of aluminum alloy bumper materials is shown in Table 2.4-1. As shown, 
the minimum weight alloy is 201 1-T8. Note the wide variation in CMC weights for different 
aluminum alloy bumper materials. 

Figure 2.4-1 shows the distribution of optimal bumper and wall thicknesses by hyper- 
velocity impact region for the 201 1-T8 aluminum bumper material and 2219-T87 aluminum 
wall material. Note that the optimal bumper thickness is most heavily influenced by the projectile 
melt/vaporization region, while the optimal wall thickness is most heavily influenced by the 
projectile shatter region. 
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Figure 2.4-2 shows the percentage area under the velocity probability distribution for the 
201 1-T8 aluminum bumper. Nearly 2/3 of the likelihood of impacts is above 8 km/sec, where 
testing is not generally attainable. 

Figure 2.4-3 shows the optimal 2011-T8 bumper thickness as a function of projectile 
diameter. This relationship is quite linear. Shown in Figure 2.4-4 is the optimal 2219-T87 wall 
thickness as a function of projectile diameter. This relationship is slightly convex. Figure 2.4-5 
gives the minimum module weight (normalized to the baseline case) as a function of projectile 
diameter for the 201 1-T8 bumper case. 

Table 2.4-1. Comparison of Aluminum Alloy Bumper Materials 


ALUMINUM 
ALLOY BUMPER 
TYPE 

OPTIMAL 

BUMPER 

THICKNESS 

(CM) 

OPTIMAL WALL 
THICKNESS (CM) 

MINIMUM CMC 
WEIGHT (KG) 

2219-T87 

0.46 

0.65 

5715 

1100-H18 

0.50 

0.64 

5839 

2011-T8 

0.46 

0.64 

5665 

2014-T6 

0.44 

0.71 

5910 

2024-T81 

0.44 

0.72 

5929 

5005-H18 

0.49 

0.64 

5760 

5050-H38 

0.49 

0.64 

5768 

5052-H38 

0.49 

0.65 

5748 

5056-H38 

0.49 

0.66 

5762 

5083-0 

0.53 

0.65 

5978 

5086-0 

0.55 

0.65 

6059 

5154-H38 

0.49 

0.65 

5769 

5357-H38 

0.48 

0.64 

5737 

5456-0 

0.52 

0.65 

5942 

6061-T6 

0.48 

0.64 

5695 

6063-T6 

0.48 

0.64 

5737 

6101-T6 

0.49 

0.64 

5760 

6151-T6 

0.48 

0.65 

5719 

7075-T6 ' 

0.43 

0.71 

5858 
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OPT. BUMPER THICK. 
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58 . 1 % 
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OPT. WALL THICK. 


Figure 2.4-1. Optimal Design Distribution By Impact Region (2011-T8 Bumper) 



Figure 2.4-2. Impact Velocity Distribution By Region (2011-T8 Bumper) 
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Projectile Diameter (cm) 

Figure 2.4-5. Minimum Module Weight vs Projectile Diameter (2011-T8 Bumper) 
2.5 Bumper /Wall Separation 

Figure 2.5-1 shows the decreasing relationship between minimum CMC weight and 
bumper/waU separation. The CMC weight shown is normalized to the baseline minimum weight 
of 5665 kg given in Table 2.4-1. Note that increasing the bumper/wall separation from 10 to 
15 cm results in a 25% decrease in CMC weight. The optimal bumper/wall separation of roughly 
200-250 cm which minimizes the normalized minimum CMC weight is shown in Figure 2.5-2. 
Finally, the optimal bumper and wall thicknesses as functions of bumper/wall separation are 
given in Figure 2.5-3. This depicts a fairly constant optimal ratio between bumper and wall 
thickness. 
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Figure 2.5-2. Minimum Module Weight vs Bumper/Wall Separation (2011-T8 Bumper) 

Bumper 

Wftl 


0 5 10 15 20 25 30 

BumpenWaJI Separation (cm) 

Figure 2.5-3. Optimal CMC Thicknesses vs Bumper/Wall Separation 

(2011-T8 Bumper) 
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2.6 Space Station Altitude 

Figures 2.6-1 and 2.6-2 show the relationships between Space Station altitude and pro- 
jectile diameter and minimum CMC weight, respectively. Note the high sensitivity of design 
weight to altitude between 200 and 1000 km. The optimal bumper and wall thicknesses as 
functions of Space Station altitude are given in Figure 2.6-3. 



Figure 2.6-1. Projectile Diameter vs Space Station Altitude 



An Employee-Owned Company 



— 



— 


2.5 

— 


2.25 


£ 

2 


? 



5 

1.75 


§ 

o 

1.5 


C 

I 

1.25 


S 

1 


N 

— 

1 

0.75 

- 


0.5 



0.25 

— 


0 

C 


im mu i 





44 


2.7 Risk Considerations 

Particle Velocity 

Figure 2.7-1 shows the normalized debris velocity probability distribution for the Space 
Station at 28.5 degrees inclination. Note the wide distribution of potential impact velocities 
from 0 to roughly 15 km/sec. Recall, also, the widely differing structural responses, and thus, 
optimal designs, over this velocity range. Figure 2.7-2 shows the cumulative normalized 
velocity probability distribution for the Space Station. 



Velocity (km/sec) 


Figure 2.7-1. Normalized Velocity Probability Distribution For 28 J Degrees 

Inclination 
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Figure 2.7-2. Cumulative Normalized Velocity Probability Distribution For 28.5 


Degrees Inclination 
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Particle Impact Angle 

The relationship between particle impact angle and velocity as prescribed by [17] is shown 
in Figure 2.7-3. Uncertainty bands are included as dashed lines. Figure 2.7-4 shows the 
normalized angular probability distribution for the Space Station. Again, the optimal protective 
structures designs vary greatly over this range. Figure 2.7-5 shows the cumulative normalized 
angular probability distribution for the Space Station. 



Figure 2.7-3. Projectile Impact Angle From Normal of Surface Oriented Parallel to 

CMC Velocity Vector vs Impact Velocity 
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Impact Angle From Normal (Degrees) 


Figure 2.7-5. Cumulative Normalized Angular Probability Distribution For 28.5 
Degrees Inclination For A Surface Oriented Parallel to CMC Velocity Vector 
Particle Arrival Time 

The panicle arrival times are generally assumed to be Poisson distributed. Thus, the 
particle interarrival times are exponentially distributed. However, the mean times of arrivals 
change over time, and therefore, particle arrival times follow a nonstationary Poisson process. 
The obvious risk associated with the particle arrival times is not knowing when impacts will 
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occur. Sensor data could reduce this risk. 

Mission Risk 

Figures 2.7-6 and 2.7-7 show the relationships between total CMC mission risk and 
projectile diameter and minimum CMC weight, respectively. The weight shown is normalized 
to the baselined weight of 5665 kg. CMC mission risk is defined as one minus the total CMC 
probability of no penetration. Note that an increase form 0.03 to 0.05 in mission risk results in 
a 30% protective structures design weight reduction. The optimal bumper and wall thicknesses 
as functions of mission risk are given in Figure 2.7-8. 
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Figure 2.7-8. Optimal CMC Thicknesses vs Mission Risk (2011-T8 Bumper) 
Mission Duration 

Figures 2.7-9 and 2.7-10 show the relationships between Space Station beginning year 
of operation and projectile diameter and minimum CMC weight, respectively. Note the convex 
shape between 1995 and 2000 followed by a concave representation through 2005. This is due 
to a benign solar flux effect in the latter years. A schedule delay of 5 years results in a 50% 
increase in protective structures design weight. The optimal bumper and wall thicknesses as 
functions of first year of operation are given in Figure 2.7-11. 
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Figure 2.7-10. Minimum CMC Weight vs First Year of Operation (Aluminum Alloys) 
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First Year of Space Station Operation 

Figure 2.7-11. Optimal CMC Thicknesses vs First Year of Operation 


Bumper 

Wall 


(2011-T8 Bumper) 


Figures 2.7-12 and 2.7-13 show the relationships between Space Station mission duration 
and projectile diameter and minimum CMC weight, respectively. These trades are for constant 
beginning years of operation of 1995. Note the shape reversal occurring at about 15 years. This 
is due to a solar flux effect for that particular period. A 10 year increase in mission duration 
more than doubles protective structures design weight. The optimal bumper and wall thicknesses 
as functions of Space Station mission duration are given in Figure 2.7-14. 
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Figures 2.8-1 and 2.8-2 show the relationships between space debris growth rate and 
projectile diameter and minimum CMC weight, respectively. Note that the design implications 
are more severe than that indicated by the growth in projectile diameter. This is due to the fact 
that the structural response of the protective structures is a nonlinear function of projectile 
diameter growth. Additionally, note that an increase in space debris growth rate from 5% to 
8 % results in a 50% increase in minimum protective structures design weight. The optimal 
bumper and wall thicknesses as functions of space debris growth rate are given in Figure 2.8-3. 
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Note that the optimal ratio between bumper and wall is fairiy constant up to about 6% debris 
growth rate, and then decreases as the wall thickness becomes a greater influence on protective 
structures design. 



0 2 4 6 8 10 12 14 

Space Debris Growth Rate (%) 

Figure 2.8-1. Projectile Diameter vs Space Debris Growth Rate 
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Figure 2.8-2. Minimum CMC Weight vs Debris Growth Rate (Aluminum Alloys) 
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Figure 2.8-3. Optimal CMC Thicknesses vs Debris Growth Rate (2011-T8 Bumper) 

Particle Shape/Densitv 

The distribution of particle shapes for space debris in orbit is unknown. The potential 
variation in protective structures design effectiveness due to changes in particle shapes has been 
shown by hydrocode and impact test data to be relatively large. 

The particle density is generally unknown as well. It is modelled as a decreasing function 
of projectile diameter as shown in Figure 2.8-4. 
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Figure 2.8-4. Space Debris Particle Density vs Diameter 


Uncertainties in Risk Parameters 

Although distributions exist for the risk parameters, these are subject to uncertainties in 
their accuracy and development. For instance, the distribution of projectile velocities is subject 
to uncertainties. Uncertainties in mission risk may be measured by establishing confidence 


intervals about the expected mission risk. 
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2.9 Second Order Parametric Analyses 

This section includes numerous design trade parametrics to aid the designer in 
decision-making and design consequences of environment-related issues. The four independent 
variables shown are bumper/wall separation, space debris growth rate, CMC mission duration, 
and CMC mission risk. 

Bumper/Wall Separation 

Figures 2.9-1 through 2.9-3 show the effects of bumper/wall separation on minimum 
CMC weight for various space debris growth rates, CMC mission durations, and CMC mission 
risks, respectively. Note, for instance, that the protective structures designer can maintain 
equivalent weight if the space debris growth rate is actually 7% by increasing the bumper/wall 
separation from 10 to 15 cm. 

Space Debris Grow th Rate 

Figures 2.9-4 through 2.9-6 show the effects of space debris growth rate on minimum 
CMC weight for various bumper/wall separations, CMC mission durations, and CMC mission 
risks, respectively. Note, for instance, that the protective structures designer can maintain 
equivalent weight if the space debris growth rate is actually 9 % by increasing the mission risk 
from 3% to 5 %. 
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Figure 2.9-1. Minimum Core Module Weight vs Bumper/Wall Separation for Various 
Space Debris Growth Rates (2011-T8 Aluminum) 



-- 10 Year Duration 

1 5 Year Duration 

20 Year Deration 

— — 5 Year Duration 


Figure 2.9-2. Minimum Core Module Weight vs Bumper/Wall Separation for Various 

CMC Durations (2011-T8 Aluminum) 
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Figure 2.9-3. Minimum Core Module Weight vs Bumper/Wall Separation for Various 

CMC Mission Risks (2011-T8 Aluminum) 
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Figure 2.9-4. Minimum Core Module Weight vs Space Debris Growth Rate for Various 
Bumper/Wall Separations (2011-T8 Aluminum) 
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Figure 2.9-5. Minimum Core Module Weight vs Space Debris Growth Rate for Various 
CMC Mission Durations (2011-T8 Aluminum) 
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Figure 2.9-6. Minimum Core Module Weight vs Space Debris Growth Rate for Various 

CMC Mission Risks (2011-T8 Aluminum) 
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CMC Mission Duration 

Figures 2.9-7 through 2.9-9 show the effects of CMC mission duration on minimum CMC 
weight for various bumper/wall separations, space debris growth rates, and CMC mission risks, 
respectively. Note, for instance, that if mission duration increases from 10 to 15 years, the 
protective structures designer can maintain equivalent weight by increasing the bumper/wall 
separation from 10 to 20 cm. 

CMC Mission Risk 

Figures 2.9-10 through 2.9-12 show the effects of CMC mission risk on minimum CMC 
weight for various bumper/wall separations, space debris growth rates, and CMC mission 
durations, respectively. Note, for instance, that if mission risk increases from 3% to 10%, the 
protective structures designer can afford to reduce the bumper/wall separation from 10 to 5 cm 
while maintaining weight. 


•s 



10 cm Separation 

15 cm Separation 

20 cm Separation 

5 cm Separation 


Figure 2.9-7. Minimum Core Module Weight vs CMC Mission Duration for Various 
CMC Bumper/Wall Separations (2011-T8 Aluminum) 
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5% Growth Rate 

7% Growth Rate 

9% Growth Rate 

3% Growth Rate 


Figure 2.9-8. Minimum Core Module Weight vs CMC Mission Duration for Various 
Space Debris Growth Rates (2011-T8 Aluminum) 



3% Mission Risk 

5% Mission Risk 

10% Mission Risk 

1% Mission Risk 


Figure 2.9-9. Minimum Core Module Weight vs CMC Mission Duration for Various 

CMC Mission Risks (2011-T8 Aluminum) 
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10 cm Separation 

15 cm Separation 

20 cm Separation 

5 cm Separation 


Figure 2.9-10. Minimum Core Module Weight vs CMC Mission Risk for Various 
Bumper/Wall Separations (2011-T8 Aluminum) 



— 5% Growth Rate 

— 7% Growth Rate 
9% Growth Rate 

— 3% Growth Rate 


Figure 2.9-11. Minimum Core Module Weight vs CMC Mission Risk for Various Space 

Debris Growth Rates (2011-T8 Aluminum) 
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— — 10 Year Duration 

1 5 Year Duration 

20 Year Duration 

— — 5 Year Duration 


0 0.05 0.1 0.15 0.2 

Mission Risk 

Figure 2.9-12. Minimum Core Module Weight vs CMC M issi on Risk for Various CMC 

Mission Durations (2011-T8 Aluminum) 
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2.10 Conclusions and Recommendations For Section 2 
Conclusions 


I, Global analytic nonlinear design optimization can be performed for the projectile 
np^a^riaifioft region (Wilkinson), for ubnn$;m^^ 

(Burch), and for the Nysnuthpi^^ 

X For the predictors investigated, the optimal rap of bumper mass per unit area to total 
mass per unit area may vary with mission, environment, projectile mass, and velocity regime. 

3. Them is a large incentive for increasing the bumper/wafl separation from 10 to 15 cm for 
all predictors investigated. ' ; , -v " 

4. M predictors fepct increasing design spsitivi^ and decreasing 

design sensitivity to bumper/wall separation. : 

5. The WilldnSOn and Ny smith predictors reflect increpng design sensitivity to projectile 
velocity, while the Burch relationship is decreasing. \ . . ,.r. : 

6, For the combined predictors, 20U~T$ is the preferable aluminum alloy bumper choice 
for the baseline parameters. V , , , : .. .. \ : ' 

7; For the combined predators, increasing the bumpcr/wall separation from 10 to 15 cm 
reduces the minimum module weight by 25%. . ” ’ 

8, Minimum CMC weight is very sensitive to spacedebns growth rate above 7% and Space 
Station altitude below. 1000 km fpr the Combined predictors. ' • •III 

9. CMC protective structures design depends greatly on mission duration for the combined 
predictors, J HP * ' 

10. For the combined predictors, increasing the CMC mission risk from 3% to 5% reduces 
die minimum module weight by about 30%. 

Recommendations 

1. Alternate metallic bumpermaterials should be investigated. V " . ' . 

X Uncertainty analyses should be performed relative to the space debris environment 
parameters. 

3. A combined meteoroid/space debris optimization algorithm should be implemented. 

4, Advanced materials should be investigated, '> - • • • : . 
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3 PROTECTIVE STRUCTURES DESIGN OPTIMIZATION CODE 

(PSDOC) OVERVIEW 

PSDOC (Protective Structures Design Optimization Code) was developed under NASA- 
MSFC Contract NAS8-37378 "Optimization Techniques Applied to Passive Measures for In-Orbit 
Spacecraft Survivability". The purpose of PSDOC is to provide a user-friendly PC environment 
for a number of design and analytical tools including IM PACT 10V, developed by SAIC. Specific 
analysis areas for spacecraft protective structures design optimization include selection of envi- 
ronment, spacecraft characteristics and mission, and hy per ve locity imp act predict or models. The 
significant features of PSDOC are a menu-driven scenario and input capability, post-processing 
features, and file management system. ~~ 

The application of S AIC’s Flexible Model - Graphical User Interface to PSDOC was but one 
utilization of this software. The graphical user interface (GUI) environment used for PSDOC was 
developed for assisting technical personnel in gaining access to co mp uter based models without a 
thorough knowledge of the code itself. Other applications are easily fitted to existing models by 
SAIC engineers and software scientists. Attachments of this GUI software to existing models or 
"Retrofitting" allows for newer coding techniques and hardware technology advancements to be 
immediately available to older, validated models without affecting the code’s reliability. Once this 
initial connection has been made and checked with the original version, additional input and output 
alterations to the model are often desired and can be handled by SAIC staff under the direction of 
our customers. 

The PSDOC environment (retrofit to IMPACT 10) was developed in coordination with 
Sherman Avans and Jennifer Robinson of NAS A-MSFC and Robert Mog, Andy Laidig, and Kevin 
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Leonard of SAIC. The PS DOC user’s manual delivered to NASA-MSFC in Aug. 1990 presents 
an overview of the windowing techniques and operating instructions for use of the PSDOC envi- 
ronment. This manual contains all the necessary information for efficient use of PSDOC. 
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4 MONTE CARLO SIMULATION ANALYSIS TOOL 

4.1 Monte Carlo Simulation Purpose 

The purpose of this simulation is to provide a statistical tool to address and quantify 
protective structures design risks, uncertainties, and options, and to address system-level issues 
relevant to designer decision-making and possible implications. The system of initial interest 
is the structural configuration of WP01, including the Core Module Configuration. "Grow- to" 
systems include module internal configurations and external structures (trusses, solar arrays, 
etc.) as specified in the redesign. 

Initial investigations of interest include statistical analyses of primary impacts, penetra- 
tions, and vulnerable areas. "Grow-to" investigations include interior effects, secondary ricochet 
effects, and SSF element interrelations. 

Risk considerations include environment particle velocity, impact angle, and component 
probability of impact. Uncertainty considerations include SSF IOC/FOC, particle diameter, 
mass-density, shape, and uncertainties in particle velocity and impact angle distributions. 

4.2 Monte Carlo Simulation Development Approach 

The tool development approach is to define the current SSF mission parameters and design 
configuration, and interpret the geometry mathematically using FASTGEN. The mission 
parameters drive requirements specification, including environment definitions. These con- 
siderations, combined with appropriate random number modules and the FASTGEN results, 
produce the necessary shotline time histories and intersecting body calculations. Survivability 
assessments follow and employ deterministic models for hypervelocity penetration prediction. 
Statistical assessments follow to supply answers to the questions of interest. 
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The top-level version of the Monte Carlo simulation tool will be executed on IBM- 
compatible PC’s. The current version of this tool runs on a VAX. It is anticipated that the 
detailed version of this tool will operate on a CRAY. I/O requirements are discussed in Section 
4 . 4 . 

Verification and validation of the tool will be performed using a combination of PSDOC 
and BUMPER. If the program execution times are considerable, a design of experiments 
approach will be used to specify production run matrices. 

4.3 Particle Time-Arrival Process for Monte Carlo Simulation Development 

Several algorithms have been developed for the particle time-arrival process. The standard 
assumption in this area is that arrival times are Poisson distributed. This means that the inter- 
arrival times are exponentially distributed, and sorting of arrival times is not required. Mean 
data is derived from the environment flux and appropriate spacecraft areas. This algorithm 
leads to a terminating simulation defined by the mission profile. 

Realistically, however, the meteoroid and debris environments are both nonstationaiy 
Poisson processes, at best, since the mean arrival rates vary in time over the mission profile. 
An approximation algorithm has been developed which alters the mean arrival rate to represent 
the time period under consideration. However, this algorithm is not exact, since a period of 
high arrival rates could be neglected using a low arrival rate corresponding to the previous 
period, or vice versa. Thus, a more exact (continuous) algorithm should be developed. The 
approximating algorithm for the space debris environment is given as: 
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If independent mean and variance data for arrival rates are available, a uniform arrival 
process may be used as an alternative to Poisson arrivals. To compare this approach with the 
Poisson process, the variance may be set equal to the square of the mean. An algorithm has 
been developed for independent mean and variance data. 
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Augmentation/repair times may be modelled using a number of distributions, if this 
modelling is of interest If mean data only for time to repair is available, an exponential service 
model may be used. If independent mean and variance data are available, the gamma, weibull, 
lognormal, or beta distributions may be appropriate. 

4.4 Simulation Status 

To date, the following items have been completed: 


A. Enveloping Geometries Established For. 

1. Sphere: Enter Radius " ' ' 

1 Cylinder: Enter Radius, Length • ... . ' ' _ * M 

3 . Box: Enter length. Width, Height - " ... 

B. Konstation^ry Poisson Arrival Process Algorithm Few Space Debris ; 

L First Random Variate Establishes Point On Cumulative Flux-Diameter Curye At 
Current Mission Time . .. | 

' 2, Absolute Flux Is Inverted To dive Mean Interarrival Time ' .4 • . 

3, Second Random Variate Establishes Time Between Arrivals Using Exponential 

Distribution ' . 

4, Cumulative Flux-Diameter Curve Is Updated For New Mission Time 

C. Impact Characteristics For Space Debris ... 

L Impact Velocity 4 ... .. 

2. Impact Angle 

3. Particle Density/Mass 

D. Look-Up Tables , \ 

L Solar Flux (Monthly) 

2. Inclination Factor . .. ' 

3, Flux-Diameter Curves 4 ... . 

EL Impact History Data Including Event Time, Diameter, Density, Mass, Velocity, 
Angle ' ' - r . ‘ ' ' 

F. Fixed Time Data Including Absolute Flux, Normalized Flux, Cumulative Nor- 
malized Flux Distributions As Functions of Diameter 
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The following items remain to be completed: 


t Impact Location/Orientation " - ' s x _ v , \ ' 

IL Meteoroid Environment . J, ; L fp? ~ 


ju* vjvOM4Cuy/ 0 noumq <iuqu v -, ^ §>? 

IV. Graphical Outputs ^> v ' \ 4 ^ : - . 

, A. Time/Interarrival Time Histories and Distributions 
; B. Particle Diameter/Mass, Velocity Distributions 
C Particle Location Display/Contours ' : , * 

D. Statistics Modules ;:; ‘ ,1/' ' .'1 

~ Plots, e.g. Diameter v|. Velocity vs. Time : f ' ; . - . 

• :• • • ; “ 


79 


5 DEVELOPMENT OF ADVANCED SHIELDING CONCEPTS 

5.1 Introduction 

The development of advanced shielding concepts presented in this section is a preliminary 
theoretical modification of the Wilkinson and ballistic PEN4 predictors to multiple bumper 
situations. The intent is to perform this preliminary analysis, and then correlate the results with 
existing test data to improve the models. 

5.2 Extension to Multiple Bumpers for Wilkinson Predictor 

A number of different approaches have been attempted to modify the Wilkinson predictor 
mathematically for multiple bumper systems. The one successful approach (physically) is given 
as follows: 

1. Modify the Wilkinson form in a product sense as: 

0.364D 3 p V cos(0) , Dp, _ 

'.= — for rn— £1 > I 981 


mu? |p. 


n p^ 

i-l 


0.364D p^V cos(0) „ 
t„ = — t — : — T7T-. r— for 




_£p£_ 

Jl-l 

n p* 


>i. 


[99] 


If our goal is to minimize system mass per unit area subject to the total separation between first 
bumper and last wall equal to some desired value, we may write this as 


-i O. 364 D 4 o*Vcos 0 


min W = X m, + / , — v 

7T- T^j 

[ 100 ] 

‘■ l dru? 

nm. 


V * 1 ) 

v-» ) 


n - 1 



S J* X = SjoT 
i * 1 


[ 101 ] 

where m, = p,/, 


[ 102 ] 
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St„ t is the total separation between the first bumper and the wall, and n-1 is the total number of 
bumpers (n is the total number of plates). 

Under condition [99], the dual Geometric Programming objective function is given by 

j, 


max v(5) = n (1/5, )*' (K75,) 8 ' 

i-i 




O.364D 4 p*Vcos(0) 


K = 


[103] 

[104] 


15 , = 1 

i ■ t 


5, -8, = 0, i » 1,2,...,«-1 

-25, + 5’ = 0, ; = 1,2,..., n-1 

m-.EA 

J m l 


[105] 

[106] 

[107] 

[108] 


Note that the degree of difficulty is 0, with 2n-2 independent variables corresponding to the n- 1 
bumper areal densities and the n- 1 separations. 

Equations [106] and [107] together imply 


i = 1,2 — 1 

[109] 

j 1,2 1 

[110] 


The minimum weight and globally areal densities are given by 


W n = n 


' 0.364D 4 pJV cos(0)‘ 


K 

l S ToT ) 


m *# = 


O.364D 4 pJVcos(0) 


-l Vn 


2* -2 


n-1 

J V. ^ToT 


i = l,2,...,n 


[ 111 ] 


[ 112 ] 
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5.3 Results 

Several results using the development of Section 5.2 are given in this section. The baseline 
assumptions are a particle density of 2.8 gm/cm 3 , velocity of 9 km/sec, diameter of 1 cm, 
impacting normally into a configuration with a total bumper/wall separation of 10 cm. 

Figures 53-1 and 53-2 show how the optimal protective structures design configuration 
varies with number of bumpers for projectile diameters of 1 and 3 cm, respectively. Note that 
for a 1 cm particle diameter, the optimal number of bumpers is 2, while for 3 cm, it is 3 bumpers. 
Also, note the significant penalty for choosing the wrong number of bumpers in these cases, as 
well as the lack of symmetry of these penalties about the optimal number of bumpers. 
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Figure 5.3-3 shows the optimal protective structures design configuration including 
optimal number of bumpers as a function of particle diameter. Increasing particle diameter 
results in an increasing optimal number of bumpers to defeat the particle. Note the optimal 
transition regions between 1 and 2 bumpers (corresponding to particle diameters between 0.75 
and 1 cm) and 2 and 3 bumpers (corresponding to particle diameters between 2.25 and 2.5 cm). 
Also, note the very linear minimum system areal density, showing the stabilizing effect of 
increasing the number of bumpers in the configuration. 

Figure 5-3-4 shows the optimal protective structures design configuration including 
optimal number of bumpers as a function of particle velocity. The most striking feature of this 
trade is the relative insensitivity to velocity for a dual bumper system. 

Figure 5-3-5 shows the optimal protective structures design configuration as a function 
of total bumper/wall separation. As in previous studies, there is a large weight incentive for 
increasing the total separation. Furthermore, increased separation allows for more bumpers to 
disrupt the incoming particle. 

Figure 5-3-6 is a replica of Figure 5.3-5, except that the optimal individual separations 
are included. 
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Parade Density - 2.8 gm/cm3, Normal Impact 
Partlcki Velocity • 9 krrVsec, Diameter - 1 cm 


Figure 5.3-5. Optimal Protective Structures Design Values vs. Total Bumper/Wall 

Separation (Modified Wilkinson) 
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Pvtide Density - 2.8 gm/cm3, Normal Impact 
Particle Velocity » 9 knVsec, Diameter « 1 cm 


Figure 5.3-6. Optimal Protective Structures Design Values vs. Total Bumper/Wali 

Separation (Modified Wilkinson) 

5.4 Extension to Multiple Bumpers for Ballistic PEN4 Predictor 

The multiple bumper recursion equations are given by: 


7 

V / = 4100, ^ £0.4 

(j \ 0.21 J 

V /= 4986 UJ • 1? >0A 




The first bumper is penetrated if 


0.67} 


( 0.28 lDp. Y' 3 

Hr] cos ' e > 


\ 1/0.31 

' 25 ,/ 


[114] 

[115] 

[116] 
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V>V. 




V- 1 


[117] 


The residual velocity (from the first bumper) is 

1 .33 V 2 /? 3 p p - (85 y ,T lg -° °°” l25V ^ cos(8) 


Vk« 




1 .33/?^p, + RpTipJ cos(0) 


The second bumper is penetrated if 


V Ky > V X i.l 


The residual velocity (from the second bumper) is 


v h m 


1.33 Vi 


[ r 2 pPp~{ 


8 SyT 2 e 


-0.0003125V, \ 1 1/2 


'1 


cos(0) 


\33R;p„ +/J,r 2 p 2 /cos(0) 


The third bumper is penetrated if 


V D >V. 


*1 ' 50 /. 3 


[118] 


[119] 


[ 120 ] 


[ 121 ] 


The residual velocity (from the (n-l)st bumper) is 


V R = 


, , ( -0.0003125V, 'N 1 ,/2 

l-33Vy ,*,P, -|cos(8) 


The nth bumper is penetrated if 


1 33Rlp p + R p T n -ip H -x/ cos(0) 


V R > Vjo 


[ 122 ] 


[123] 


Given 6061-T6 aluminum bumper materials (yield strength of 35 ksi, density of 2.71 
gm/cm 3 , total thickness of 0.16 cm), 2219-T87 aluminum wall (yield strength of 51 ksi, density 
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of 2.81 gm/cm 3 , thickness of 0.3175 cm), a projectile density of 2.81 gm/cm 3 , and a projectile 
impact angle of 0 degrees (normal). Figure 5.4-1 shows the ballistic limit curves for single, 
double, and triple bumper configurations. Note the relatively minor sensitivity to number of 
bumpers over this limited range. 



Single Bumper Double Bumper Triple Bumper 

(Two Plate) (Three Plate) (Four Plate) 

— « — — ■ — — * — 

6061 -T6 Aluminum Bumper 
2219-T87 Aluminum Wall 

2.61 gnrV'cubic cm Projectile, Normal Impact 

Figure 5.4-1 Critical Diameter vs. Projectile Velocity for Multiple Bumper Systems 

Using Ballistic PEN4 Recursion 
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5.5 Advanced Shielding Concepts Status 

To date, these multibumper concepts have been shown for a theoretical modification of 
the Wilkinson predictor, as well as for the ballistic PEN4 predictor. It is recommended that 
these concepts be extended to the Burch predictor, and that the Wilkinson extension be correlated 
with hydrocode data and the Burch extension with impact test data. 
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6 DISCRETE PROTECTIVE STRUCTURES DESIGN OPTIMIZATION 
6.1 Introduction 
Background 

Within the field of nonlinear programming lies a technique called geometric programming. 
Geometric programming is a purely algebraic method that provides global, and often analytic, 
solutions to certain problems previously discussed in Section 2.3. These problems are called 
posynomial programs, and they are generally nonconvex, nonlinearly constrained formulations. 
The field of geometric programming has been extended to programs which are not 
posynomials; 43,44,117 ' 122 however, the global features of the solution are not retained in this 
extension. Thus, the term posynomial programming, sometimes called prototype geometric 
programming, refers only to programs composed entirely of posynomials. 

In general, discrete nonlinear optimization techniques are even less capable than con- 
tinuous ones of providing global and analytic solutions. 123 In particular, many current discrete 
nonlinear techniques employ branch and bound derivatives, which generally do not result in 
global optimization properties, except for convex programming problems. Discrete posynomial 
problems which can be transformed to prototype geometric programs, on the other hand, result 
in global optimization upon transfer to the dual. The general transformation can then be applied 
to engineering design problems with independent variables restricted to standard or discrete 
availabilities. 

Subtask Goal 

This subtask addressing the development and application of discrete nonlinear optimi- 


zation techniques is not required in the Statement of Work, but is a natural extension of the 
traditional continuous optimization problem. A full treatment of this subtask is given in Discrete 
Posynomial Programming With Applications To Spacecraft Protective Structures Design 
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Optimization , by R. A. Mog. The goal of this subtask is to develop a theory for solving nonlinear 
programming problems that may be stated in standard posynomial form under the guidelines 
of prototype geometric programming, but with discrete constraints on the primal independent 
variables. The main development thrust is in the direction of dual program solution methods, 
although primal solution techniques are also developed. Dual method solution approaches will 
depend on problem degree of difficulty, but for problems with nontrivial degrees of difficulty, 
partial invariance and direct search techniques are investigated for their utility. Because sig- 
nomial (polynomial with undetermined coefficients and real exponents) programming methods 
do not result in global optimization, extensions of discrete techniques to signomial and reversed 
inequality constraint problems are considered secondary to this effort. 

Another goal of this subtask is to demonstrate applications of the developed discrete 
posynomial programming methodologies. These applications include challenges in the field of 
spacecraft protective structures design optimization and emphasize missions that are relevant 
to Space Station Freedom and space debris/meteoroid environments. Specific hypervelocity 
impact predictor models include those of Nysmith, Wilkinson, and Burch. 

Subtask Approach 

After a brief review of posynomial programming in Section 6.2, two primal methods for 
solving the discrete posynomial program are introduced in Section 6.3. The methods are 
numerical in nature and easy to apply to practical problems. However, no global or analytic 
information is guaranteed in their application. Therefore, in Section 6.4, a dual method is 
developed which provides the global optimal solution to the discrete posynomial program. 
Finally, three case studies which illustrate the capabilities of the primal and dual methods in 
this field are presented in Section 6.5. 
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6.2 Review of Posynomial Programming 

In a search performed at the Redstone Scientific Information Center (RSIC) to determine 
documents with the keywords "Geometric Programming '' m either the title or abstract, a total 
of 92 listings were found. Of the 92 listings, approximately 34 were theoretical and 58 applied. 
Most of the theoretical listings dealt with algorithmic improvements, code comparisons, tutorial 
papers explaining the method, and theses on specific areas of geometric programming devel- 
opments. * Of the 58 applied listings, almost all involved structural design appli- 
cations. u,17 ’^' 58 ’ 70 ' 79 ' 92 ’ 96,107 ' U2,IM,133 ' l4i ' 151 ' 152 ' 155 Other applied areas included economic,* 7 
communications, 69 and traffic flow problems. 59 Perhaps most surprising is that 27 of the 92 
listings were written after 1980. Since geometric programming was formalized in 1967, this 
points to a possible resurgence in the method’s use. 

The most interesting conclusion from the article survey is that the relatively large number 
of application papers conflicts with the dismissals of many textbook authors concerning the 
utility of geometric programming. With this many listings, it is clear that some scientists are 
finding great uses for GP. 

Based on the relatively large number of applied geometric programming listings in the 
article survey, it is apparent that GP possesses a fairly high utility, particularly in the area of 
structural design. Because GP is the only nonlinear programming (NLP) technique which offers 

* . . * : ...i. .... 

the guarantee of a globally optimal solution for certain nonconvex problems, it should be 
considered more widely in practice. Additionally, for zero degree of difficulty problems, GP 
can provide an analytic optimal solution for the objective function and independent variables. 
This attribute provides greater insight for the system designer than that obtainable by other NLP 
techniques. Finally, the values of the dual variables may provide very crucial design information 
alone in terms of the physical parameters of the problem at hand. 
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6.3 Discrete Posynomial Programming Using Primal Methods 

Introduction 

As explained earlier, geometric programming includes both posynomial and signomial 
programs and dual and primal approaches. In this section, a primal method for solving discrete 
posynomial programs is developed. The technique employed is an exterior penalty function 
method 10,53 supported by two alternate search techniques: a random/exhaustive search 53 and a 
Hooke and Jeeves pattern search . 10,53 The primal methodology computer code is called 
POLYPRIME.FOR and is given in Appendix A. 

Penalty Function Development 

Penalty function techniques are widely used numerical optimization methods which 
convert constrained optimization problems into unconstrained ones with appropriate penalties 
for not satisfying the constraints . 10,53 Two general classes of penalty functions exist: exterior 
and barrier functions. Exterior penalty function methods generally begin with points outside 
the feasible solution space and progressively drive the solutions into the feasible region. Barrier 
function methods require feasible initial points in setting up blockades along the constraint 
surfaces. 

For the problem of solving the primal formulation of discrete posynomial programs, an 
exterior penalty function technique is chosen to relieve the analyst of the burden of specifying 
a feasible initial point. This requirement could be particularly difficult when combinations of 
continuous and discrete constraints are involved. The primal problem is specified as 
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and the additional discrete constraints 

Xj — Kjfij, j — 1, 2 , . . ., (j ^ k 

One choice for an unconstrained penalty function is 


[126] 


min<|> = X c, II x?+ £ 8^4, 

i-J >-l /-i 
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where 


5 ,= 1 , g-K,> 0 

8/ = 0 > g/—Ki< 0, / = l,2,...,p 

A, = l, y = 1,2,...,<7. 


[128] 

[129] 

[130] 


Here, it is assumed, without loss of generality, that the Xj’s may be reordered such that the first 
q of them are those requiring discrete solutions. Also, the discrete penalty term has an exponent 
of 1/2 to require a stricter measure of convergence, since 


|< 0.5 < 1 


[131] 


In POLYPRIME.FOR, the accelerating factors begin at 1.0 and progressively are multiplied by 
10 until convergence is reached, i.e. 


i 5 A l Cu n i fi-\ x f] i l/ 2 < e =o.ooi 

/=1 V-i "i-1 7 ) ; = 1 J rj 


[132] 


Note that this method handles mixed discrete problems as well as continuous and purely discrete 
problems. Furthermore, note that although we are strictly concerned with posynomial pro- 
gramming problems, this penalty function approach is equally valid for generalized polynomials 
or signomials and signomial constraints. Additionally, Type I, II, or in inequality constraints 
are valid. However, the constraints must be converted to Type I, less than or equal to constraints. 



An Employee-Owned Company 


97 


for implementation in POLYPRIME.FOR. Similarly, the constraint right hand side values, K„ 
may take any real value. This is a more generalized form than that allowed for prototype 
posynomial programming where the K,’s must be equal to 1 for all constraints. Finally, note 
that if convergence of continuous equality constraints is difficult to achieve using this formu- 
lation, it may be easily modified by adding a penalty term for equality constraints separate from 
inequality constraints. Note, however, that continuous equality constraints combined with 
discrete variable constraints could easily result in no feasible solution situations. 

Now, once the unconstrained penalty objective function is established, a method to 
minimize it must be found. Two approaches using search techniques are discussed in the next 
two sections, followed by a comparison of the methods. 

Random/Exhaustiv e Search Subtechnioug 

A random search technique 53 is analogous to throwing darts at a dartboard with no 
adaptation or learning between throws. (There do exist adaptive random search techniques, but 
these won’t be considered here.) Although random search techniques may appear unsophisti- 
cated due to their brute force nature, they are particularly useful in establishing optima of highly 
nonlinear and multimodal functions. Since discrete nonlinear optimization problems tend to 
add a degree of this type of complexity, it would appear that random search techniques would 
prove fruitful. 

The number of search points for a pure random search is given by Gottfried as 

m, = iln(l-P,) I 133 ] 


The search space is defined by specifying an interval of interest for each of the independent 
variables. The main drawback for employing a random search technique in a discrete optimi- 
zation problem using penalty functions is the severity of the convergence criteria. Unless the 
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random draw is extremely fortuitous, the convergence criteria will not be met, particularly for 
sparsely populated discrete feasible regions. For this reason, an exhaustive search option is 
automatically called in POLYFRIME.FOR when the number of discrete feasible points (as 
specified by the search space) is less than the number of search points given by n\ above. On 
the other hand, if the search region is dense with discrete points^ as compared with the pres- 
pecified number of random search points, then random search proceeds normally. 

Hooke and Jeeves Pattern Search Subtechniaue 

The Hooke and Jeeves pattern search 1043 is a more methodical unconstrained search 
technique, which requires an initial point, but no variable search intervals. The technique begins 
with exploratory moves to establish a base point. These moves are followed by pattern moves 
through successive base points. Convergence requirements are more easily met for discrete 
problems using this method. 

Comparison of Subtechniques 

The Hooke and Jeeves pattern search technique is more methodical and generally con- 
verges faster than the random/exhaustive search technique. Furthermore, it requires only an 
initial point rather than an interval of investigation. On the other hand, the Hooke and Jeeves 
method is generally fairly sensitive to the initial search point and is less likely to find global 
optima for multimodal functions. Furthermore, although the random/exhaustive search tech- 
nique may overly restrict the region of interest for a variable, this condition can easily be 
diagnosed when the optimal solution is found at an interval endpoint. 
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One excellent approach to combining the two methods’ strong points is to use the random 
search technique in solving the corresponding continuous optimization problem, and then use 
that solution as the initial search point for the Hooke and Jeeves subtechnique for the discrete 
problem. Another interesting study would be to compare the two methods from a time- 
effectiveness standpoint. 

6.4 Discrete Posynomial Programming Using Dual Methods 

Introduction 

The use of penalty function techniques with random/exhaustive and Hooke and Jeeves 
search subtechniques may provide rapid convergence for discrete posynomial (and signomial) 
programs. However, numerical instabilities may occur in the penalty function acceleration 
parameters due to ridges in the penalty function. Additionally, global optimal solutions are not 
generally achieved, particularly when exhaustive search is not performed. Finally, little ana- 
lytical information is gained for sensitivity analysis when numerical methods are applied. The 
restriction to posynomial programs does not lead to any significant advantages over signomial 
programs using the primal approach. Indeed, the method is perfectly valid for generalized 
polynomials or signomials and general constraints. 

In this section, dual methods are applied to discrete posynomial programs. The main 
advantages to these approaches are the guarantee of a global optimal solution and the analytic 
information gained during the process. The main disadvantages are that some derivation is 
required and that many nondegenerate continuous programs result in discrete programs with 
high degrees of difficulty. 
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General Development and Degree of Difficulty 


The prototype discrete primal program is given as 


min/=Xc,n^ 

1-1 Jm 1 

[ 134 ] 

g , = x c u n x‘* <. i, /= 1,2, ...,p 

i - 1 j = 1 

[ 135 ] 

Xj = rjn jt j = 1,2,..., <7 

[ 136 ] 


where c it Xj, c u , and ij are positive valued for all i, j, and 1, and n^ is a positive integer for all j. 


Note that this primal program is quite similar to that given in Section 6.3, with the only difference 
being the replacement of the K,’s with the value 1 . Additionally, in the case of the dual methods, 
the positivity restrictions are strictly adhered to. Recall that the primal method is equally valid 
for all real coefficients, general right hand side values, and Type I, ff, and in constraints. 

The first and most obvious fact to note about this problem is that the continuous optimal 
objective function value (obtained by not considering the discrete constraints) is always a lower 
bound for the discrete objective function value. This is easily seen by contradiction, since, for 
any set of nj’s established in a discrete optimal solution, the independent variables may always 
take the values n^ for any equivalent continuous problem. 

The second item to notice is that the discrete equality constraints may always be written 
as pairs of prototype posynomial constraints, i.e. 

Xj rj l: 

— ^Ia-^^1 [137] 

r j n j x j 


Based on this observation, we may establish a first theorem for discrete posynomial program- 
ming. 
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Theorem 1. Suppose {n 1 ,n 2 v»n q } are positive integers. Then, provided the dual program 
is consistent and the feasible discrete solution space is nonempty, the mixed discrete posynomial 
programming solution is globally optimal. Furthermore, the discrete dual program degree of 
difficulty is 2q greater than that for the continuous dual program. 

Proof: 

By the preceding observation, the discrete primal problem can always be formulated as 
a prototype posynomial program, which has been shown to be globally optimal (provided 
feasibility/consistency relationships hold) under the dual program 


maxv(5)= II 

t » 1 


|V 


’ [1381 


with 


t - (8a - 5„) + if 1 5>^,j = 0 h = \,2,...,q 


Z + z 

i-i ;-r 




m, 


Z 5;,a >w =0 h =q + l,q + 2,...,k 
V- 1 ) 


I5, = l 

i * 1 


[139] 

[140] 

[141] 


H,= Z5; l = 1,2, ...,p [142] 

}• * 

Thus, there are k+p+l equations in n+2q+p+m,+m 2 +...+m p unknowns. The discrete dual degree 
of difficulty is given by 

DOD =n -k- 1 +2q + £ m, [143] 

;>i 


which is 2q greater than that for the continuous dual problem given by equation [32]. Thus, the 
theorem is proved. 
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Note that this theorem is not particularly useful without some idea of the nature of the 
desired integer values nj. The following theorem is useful as an expression for the basic dual 
variables. 

Theorem 2. For the discrete dual program, the basic dual variables may be written in 
terms of the nonbasic dual variables, the nj’s, and the nondiscrete primal variables as 


5„ = c, 


n x‘« 


I 

f-l 1 




FTT 

I if . 
0-1 *)( 


) 


[144] 


form=l,2,...,n. 

Proof: 

The global optimal solution is given at the equality of the arithmetic and geometric means, 
i.e. when 

k 


c, jTI x/ = 5,v(5), i = l,2,...,n. 


[145] 


This is often referred to as the dual-to-primal transformation. Thus, we may write 


(i 


n^=n 

)!• I i-1 


5, 


V V 'V 


fUrjn,)^ ft |f 

>■1 /-I 


V 


n 

f-n 




1 = 1,2,...,/! [146] 


or 


f - A 8 *" 1 


n^= 




n F 

i-il o, 


j * 1 i = 1 


C JI 


n 


>-ii 


j 


, m = 1,2,...,/! [147] 


But 
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Xj = rjtij, j - 1 , 2 , ...,<? 


[148] 


(B m f 

UJ 


h($f J fife?] 

1 i*i\^/ y»l l m l J J 

i+m v ' 


m « l,2,...,n [149] 


or 


f 8 .r- 

bJ 


n x? 

/** + 1 


n *;■ 




m = l,2,...,n [150] 


The desired result follows by substituting 

£ 5, = 1 - 5,,, m = l,2,...,n 


[151] 


i-1 
* #« 


in the exponent. Use of this result combined with the dual linear equality constraints may help 
define further avenues of solution. 

"Posyseparable' 1 Programs a nd Partial Invariance 

In this section, the term "posyseparable" is introduced, followed by techniques which may 
simplify the solution of the dual discrete posynomial program, including partial invariance. 

Definition: A posynomial function, f, is called "posyseparable" if each of its independent 
variables may be isolated at least once, i.e., if it may be written 

* . »-* * 


f- 1 erf + 1, c, Ylx? 

1-1 l-l jm 1 


[152] 


where the x/s, c u ’s, and c ( ’s are positive valued. Note that the term posyseparable, as applied 
to posynomials, is less restrictive than the term separable (see separable programming) where 
each independent variable is completely isolated in a functional sense from the remaining 
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independent variables. Many posynomials d isplay the posyseparable property. Although this 
property is not a necessary condition for existence or uniqueness of global optimal solutions for 
discrete posynomial programs, it does allow for a straightforward dual-to-primal variable 
conversion, which is often a major drawback to using dual methods. 

Theorem 3. If the objective function in a consistent and feasible discrete posynomial 
program is posyseparable, then the basic dual variables may be written in terms of any of the q 
discrete variables as 


5*. = c w [ 


V 5 W 


5-1 


n r| n 

1-1 l O u } i - 1 




i 


Ci 

Sj 


W 




• n 

/-i 


/ 'b C V, A 

g-r-M.* 

\°V C M 


V*.; 


p 

n 


15., 


J. 1 ') 


"*( Cji 

n F 


V/r 


1- I 5 r I S L 
<•1 1 i-l w 
| i #** 


,/n = 1,2, — , A: [153] 


where 


Sj* S u — afib + X Qu,&i + X 


X 5' a. 




i i = i^,-i V 


, h = 1,2,. ..,<7 


Proof: The dual program is 


maxv(8)= FI 


fr, Y%-* 


i*i 




n 




[154] 


[155] 

^*1 ) 1 . 1 y = l 8 ;, 


Mi + 2 5 ( a A -(S 2 A-8 1 *)+ ij 

» * 1 / = 1 


X8 jflju =0 h = 1,2,. ..,<7 


V-i 


+ X 8,a * + £ X 8:,a^, =0 /t=<7 + l,<7+2, ...,i 


i-l i ■ 1 v * 1 




* »-* 

.X5„+ X5, = l 


i=l i-l 


[156] 

[157] 

[158] 
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Since f is posyseparable, 


Also, 


ft =S8^ l=l,2,...,p 

i m i 


C mt X m — 8*,V(5), m — l,2,...,k 


Xj = r j n j , j- 1,2 q 


Therefore, 


f - \ 


v 5 ~y 


k:-=n 

i * 1 




WrY' 


n 

1*1 


£i 




But 


or 


Therefore, 


a m 

r r ’ 

^<r=v( 5)=^f- 


- j — V 


'O a ~ 


f = n 

0 „ 1*1 


n 

i»i 


v 5 «y 


£i 

v 5 * y 




•(ftf: 


"i , 


/ , 

M ”, 


ai&f')), m = l,2,...,k 

j ‘ * 5 >; 


Since 


*, = V*« 



[160] 

[161] 

[162] 

[163] 

[164] 


[165] 


[166] 
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V J ' V ' 


\i 

,-1 , , ,(«| V-‘ * "I c 

.nc^A'^fi 2^ (n(^)Oi 

;-i 7 7 <-i^>-i v y-i Oj* 


m = 1,2,..., it [167] 


r n -\ fr m 

r >”> X 


fc Y* _1 * ( C .^.Vc-^ 1 5. . 

s_=c„[ f n f n j (rsJ* y • 

lO„J 1-1 l o^l «-llO, J 


[168] 


« -»[ 8;, o * A r-k *( Ca 
n 5 -— w n i 5 ; , n / 

>* ,l v*V c > J , “ t v* s J u**l8jj 


N n-k k 

AV-,!i'-.i> 


,m = \,2,...,k [169] 


m <q - 1 


[170] 


wc may rewrite 


n-k k 

5^ = 1- I 5,- 15, 

»*1 i» 1 
1 # m 


[171] 


and the theorem is proved. Combining this result with the linear dual equations allows one to 
use partial invariance to solve for basic dual variables in terms of at most one nj. Then, a series 
of differences in terms of the discrete and continuous dual objective functions may be minimized 
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v<(5)-v c (5)=^p--v c (5) [172] 

Setting this difference to zero, we may write 

- S-V t (5) = - 5„v,(S) [173] 

=«Jr (, '.' , . ) ' , - v * (8) ] 11741 

\°v J 

which is a function of n, and the dual variables only. In general, partial invariance should be 
used to solve the basic dual variables in terms of the dominant dual variables if one has that 
knowledge for the problem at hand. A good starting point is to solve for the dual variables 
corresponding to the discrete (primal) constraints in terms of the basic dual variables evaluated 
in the neighborhood of their optimal continuous solutions. 

A fair amount of derivation has indicated possible directions a dual approach may take 
in a discrete posynomial program. The concept of posyseparability has been introduced to ease 
the dual-to-primal transformation. The discrete dual program has been shown to have a degree 
of difficulty 2q greater than the corresponding continuous dual program. A number of relatively 
simple examples have been used to illustrate various facets of the dual approach. Primal and 
dual approaches have been compared. Specific numerical and computer methods used to support 
the dual approach are left for future development. 

Figures 6.4-1, 2, and 3 show the optimal discrete design values of minimum system mass 
per unit area, and optimal bumper and wall thicknesses vs projectile diameter, bumper/wall 
separation, and projectile velocity, respectively, for the Nysmith predictor. The discrete 
availability factor, r„ is 1/64 inches. Figure 6.4-1 reflects a constant meteoroid density. In 
Figure 6.4-3, the impact angle remains constant at 0 degrees (normal). 
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Figure 6.4-1. Optimal Discrete Design Value vs Projectile Diameter for Nysmith 

Predictor 
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Figure 6.4-2. Optimal Discrete Design Value vs Bumper/Wall Separation for Nysmith 
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Optimal Discrete Bumper Optimal Wall Minimum Total System 

Thickness (cm) Thickness (cm) Thickness (cm) 


Figure 6.4-3. Optimal Discrete Design Value vs Projectile Velocity for Nysmith 

Predictor 

Figure 6.4-4 shows the sensitivity of minimum system mass per unit area to bumper 
thickness availability factor, r,. The discrete and continuous objective functions are equal when 
the continuous bumper thickness is an integer multiple of the bumper thickness availability 
factor as shown in Figure 6.4-5. This occurs at numerous locations over the range considered. 
Note that when r t is small, the discrete bumper thickness is closer in value to the continuous 
bumper thickness. As r, grows, this incidence of equality naturally decreases while the deviations 
from the continuous minimum mass per unit area grow in value. Beyond the optimal continuous 
value of the bumper thickness, the objective function continues to grow indefinitely, because 
the availability factor is dominating the desired continuous solution. 
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Figure 6.4-4. Minimum Total System Mass Per Unit Area vs Bumper Thickness 
Availability Factor for the Nysmith Predictor 
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Figure 6.4*5. Optimal Bumper Thickness vs Bumper Thickness Availability Factor for 

the Nysmith Predictor 

Figures 6.4-6, 7, and 8 show the optimal discrete design values of minimum system mass 
per unit area, and optimal bumper and wall thicknesses vs projectile diameter, bumper/wall 
separation, and projectile velocity, respectively, for the Burch predictor. The bumper thickness 
availability factor is 1/64 in. Figure 6.4-6 reflects a constant projectile density as given in 
equation [141]. In Figure 6.4-8, the impact angle remains constant at 0 degrees (normal). 
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Figure 6.4-6. Optimal Discrete Design Value vs Projectile 
Diameter for Burch Predictor 
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Figure 6.4*7. Optimal Discrete Design Value vs Bumper/Wall Separation for Burch 

Predictor 
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Optimal Wall 
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Figure 6.4-8. Optimal Discrete Design Value vs Projectile 
Velocity for Burch Predictor 

Figure 6.4-9 shows the sensitivity of minimum system mass per unit area to bumper 
thickness availability factor, r t . The discrete and continuous objective functions are equal when 
the continuous bumper thickness is an integer multiple of the bumper thickness availability 
factor as shown in Figure 6.4-10. This occurs at numerous locations over the range considered. 
Note that when r t is small, the discrete bumper thickness is closer in value to the continuous 
bumper thickness. As r, grows, this incidence of equality naturally decreases while the deviations 
from the continuous minimum mass per unit area grow in value. Beyond the optimal continuous 
value of the bumper thickness, the objective function continues to grow indefinitely, because 
the availability factor is dominating the desired continuous solution. 
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Discrete Continuous 


— — 

Figure 6.4-10. Optimal Discrete and Continuous Bumper Thickness vs Bumper 
Thickness Availability Factor for the Burch Predictor 
The discrete Wilkinson predictor optimization algorithm is derived using the dual method. 
Theorem 4. The combined discrete/continuous Wilkinson algorithm is given by 


0.364D 4 p$Vcos(9) 

1 * C i ™ n 

L]S p, 

[175] 
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Otherwise, go to step 10. 


5. t,=r x n x 
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If tj = r 2 n 2 is required, then the optimal discrete values are given by 
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and return to step 5. Otherwise, 
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L 2 Pl r, L, 

Check both values and return to step 5. 

Proof: The primal objective function is given by 

W — Pi^i P 2 h ~ Pi*i +7" 

u 

and is constrained by t! = r t n,. The dual program is given by 

maxv(8)= — (w)* 1 

V 5 ‘/ V 5 */ 

s.t . 5j+6j=1 

S 1 -5 2 + 8' 1 -^ = 0 

(o^f Cl Y" 81 2*,-^ 

=> maxv(8) = -r- - — ( r i n i) 
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Now, minimizing the difference between the discrete and continuous dual objective functions 
for Wilkinson gives 
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At zero this minimum gives 
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c 1 + p,r 1 2 n 1 2 = r 1 «iW 0e 
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and the roots are identical. For 
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and the main result is proved. Note that the justification for using nearest integer is the parabolic 
fonn of the quadratic equation 
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Figures 6.4-11, 12, and 13 show the optimal discrete design values of minimum system 
mass per unit area, and optimal bumper and wall thicknesses vs projectile diameter, bumper/wall 
separation, and projectile velocity, respectively, for the Wilkinson predictor. In Figure 6.4-11, 
the projectile density varies with diameter according to equations [13] and [14]. In Figure 
6.4-13, the impact angle remains constant at 0 degrees (normal). The optimal bumper and wall 
thicknesses for the Wilkinson predictor are approximately equal due to the similarity in bumper 
and wall material densities. 



Optimal Discrete Bumper Optimal Wall Minimum System Mass 

Thickness (cm) Thickness (cm) Per Unit Area (gm/cm2) 


Figure 6.4-11. Optimal Discrete Design Value vs Projectile 
Diameter for Wilkinson Predictor 
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Figure 6.4-12. Optimal Discrete Design Value vs Bumper/Wall Separation for 

Wilkinson Predictor 
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Optimal Discrete Bumper Optimal Wall Minimum System Mass 

Thickness (cm) Thickness (cm) Per Unit Area (gm/bm2) 


Figure 6.4-13. Optimal Discrete Design Value vs Projectile 
Velocity for Wilkinson Predictor 

Figure 6.4-14 shows the sensitivity of minimum system mass per unit area to bumper 
thickness availability factor, r,. The discrete and continuous objective functions are equal when 
the continuous bumper thickness is an integer multiple of the bumper thickness availability 
factor as shown in Figure 6.4-15. This occurs at numerous locations over the range considered. 
Note that when r, is small, the discrete bumper thickness is closer in value to the continuous 
bumper thickness. As r, grows, this incidence of equality naturally decreases while the deviations 
from the continuous minimum mass per unit area grow in value. Beyond the optimal continuous 
value of the bumper thickness, the objective function continues to grow indefinitely, because 
the availability factor is dominating the desired continuous solution. 
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Figure 6.4-15. Optimal Discrete and Continuous Bumper Thickness vs Bumper 
Thickness Availability Factor for the Wilkinson Predictor 
6.5 Conclusions and Recommendations for Section 6 
Conclusions 

In conclusion, global (and sometimes analytic) optimization of discrete posynomial 
programs can be performed using dual approaches coupled with partial invariance techniques. 
However, primal methods require less "pencil and paper" effort than dual methods and are more 
easily applied to most problems. Primal methods do not generally obtain global solutions for 
the discrete posynomial program. Furthermore, the dual method may be advantageous in cases 
where the objective function may be sufficiently separable, since posyseparable programs do 
not require solutions of coupled nonlinear equations in the dual-to-primal variable transfer- 



An Employee-Owned Company 



128 


mation. For protective structures design optimization problems, global nonlinear design opti- 
mization can be performed for the Wilkinson, Burch, and Nysmith impact predictors. In these 
cases, the optimal ratio of bumper mass per unit area to total mass per unit area may vary with 
mission, environment, projectile mass, and velocity regime. Additionally, there is a large 
incentive for increasing the bumper/wall separation from 10 to 15 cm for all three predictors 
investigated. All three predictors reflect increasing design sensitivity to projectile diameter and 
decreasing design sensitivity to bumper/wall separation. However, the Wilkinson and Nysmith 
predictors reflect increasing design sensitivity to projectile velocity, while the Burch relationship 
is decreasing. 

Recommendations 

It is recommended that other primal methods be investigated, including penalty functions 
supported by derivative search methods and feasible direction developments for discrete 
posynomial programs. Additionally, computer algorithms should be implemented based on 
current dual codes and modifications to the discrete problem. The dual method should also be 
extended to signomials. In the area of spacecraft protective structures design optimization, other 
hypervelocity impact predictors should be investigated. The discrete methods developed in this 
study should also be applied to other structural design problems. Finally, alternate protective 
materials and configurations should be investigated. 
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7 HYPER VELOCITY IMPACT TEST SAMPLE DAMAGE ASSESSMENTS 

Hypervelocity impact test sample damage assessments were performed by UAH. Posynomial 
regression analysis was performed by SAIC, and is available in Discrete Posynomial Program- 
ming With Applications To Spacecraft Protective Structures Design Optimization , by R. A. 
Mog. 

The purpose of this effort is to show a posynomial regression analysis of existing hypervelocity 
impact test data, followed by the global optimization of the ensuing structural design problem 
incorporating the predictor. A posynomial (polynomial with positive coefficients and positive- 
valued independent variables, but not necessarily positive exponents) form is chosen for several 
reasons: 

1. Posynomials can be globally optimized using the nonlinear geometric programming technique. 

2. Many previously developed predictors (by Nysmith, Madden, Wilkinson, Richardson, etc.) are 
of posynomial form. 

3. Posynomial regression problems may, under certain circumstances, be solved using linear 
regression techniques, which are easier to solve and measure statistically. 

This effort focuses on the question of whether posynomial regression can be performed in a 
statistically significant manner. A secondary goal of the study is to provide global optimization of 
the design problem formulated using the derived posynomial predictor. 

The development and analysis of a posynomial hypervelocity impact predictor suitable for 
the design of protective structures for spacecraft exposed to the meteoroid and space debris environs 
is presented in the reference above. The posynomial form is first developed with a number of 
estimated parameters. This model is next transformed into a linear regression model. Regression 
analysis is performed using a least squares approach to estimate the parameters, followed by analysis 
of variance, F-tests, and correlation coefficient examination. Residual values are then plotted against 
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the predicted and variable values. Next, the model is transformed into a hypervelocity impact 
penetration predictor suitable for design. Finally, the design problem is formulated and globally 
optimized using posynomial programming. Results show that statistically significant posynomial 
impact predictors can be developed using linear regression approaches. 

The main conclusion of this effort is that it is possible to develop a statistically significant 
posynomial hypervelocity impact predictor with a fairly large number of impact tests and a fairly 
small number of predictor variables. Although greater variation can be explained by considering 
posynomials with more than one term, the ability to transform the posynomial into a form suitable 
for linear regression is lost. Furthermore, since it is generally desirable to have 10 or more data 
points per predictor variable, the increased number of term values might actually decrease the 
confidence in the predictive capability of the model as measured by the analysis of variance. 



An Employee-Owned Company 


131 


8 ANALYSIS OF PROJECTILE SHAPE EFFECTS 

SAIC developed posynomial regression techniques and combined them with posynomial 
optimization techniques for application to this area. These techniques are available for immediate 
application to the test data resulting from projectile shape effects testing. Currently, limited test 
data produces unclear results when attempts are made to correlate data from various projectile 
shapes. Results are inconclusive. Further investigation of the projectile shape effects could include 
methodologies found in sources such as "A Preliminary Investigation of Projectile Shape Effects 
In Hypervelocity Impact of a Double-Sheet Structure," by R. H. Morrison, NASA-TN-6944, August 
1972, but will remain inconclusive until further test are performed. 
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10 APPENDICES 

APPENDIX A. POLYPRIME.FOR: 

A GENERALIZED PRIMAL OPTIMIZATION TECHNIQUE USING PENALTY 

FUNCTIONS 

DIMENSION C(100)A(100,100),X(100),M(100),a(100,100) 

DIMENSION AI(100,100,100),XKL(100),XHIGH(100),XLOW(100) 

DIMENSION RJ(100),IDIS(100),XOPT(100)^EXP(100) 

DIMENSION EPS 1(100), EPS2(100) 

OPEN(UNIT= 1 0,TYPE= ’OLD ’ ,ACCESS= ’ SEQUENTIAL’) 

OPEN(UNIT=l 1,TYPE=’NEW’,ACCESS=’SEQUENTIAL’) 

♦♦♦INITIAL SEED FOR RANDOM SEARCH*** 

ISEED = 91411 

♦♦♦NUMBER OF CASES TO RUN*** 

READ(10,*) NRUNS 
DO 10 IR=1, NRUNS 

***IOPT=l FOR RANDOM SEARCH, 2 FOR HOOKE AND JEEVES*** 

READ( 10,*)IOPT 

♦♦♦NUMBER OF TERMS IN OBJECTIVE FUNCTION*** 

READ(10,*)N 

♦♦♦NUMBER OF INDEPENDENT VARIABLES*** 

READ(10,*)K 

♦♦♦NUMBER OF CONSTRAINTS*** 

READ(10,*)P 
DO 15 1=1, N 

♦♦♦COEFFICIENT FOR EACH TERM IN OBJECTIVE FUNCTION*** 

READ(10,*)C(I) 

DO 20 J=1,K 

♦♦♦EXPONENT FOR OBJ. FUNC. BY VARIABLE AND TERM*** 

READ(10,*)A(I,J) 

20 CONTINUE 
15 CONTINUE 
DO 25 L=1 ,P 

♦♦♦NUMBER OF TERMS BY CONSTRAINT NUMBER*** 

READ(10,*)M(L) 

♦♦♦RIGHT-HAND-SIDE BY CONSTRAINT NUMBER*** 

READ(10,*)XKL(L) 

DO 30 1=1, M(L) 

♦♦♦COEFFICIENT BY TERM AND CONSTRAINT NUMBER*** 

READ( 10,*)CI(I,L) 

DO 35 J=1,K 

♦♦♦EXPONENT BY TERM, VARIABLE, AND CONSTRAINT NUMBER*** 

READ( 1 0,*) AI(I,J,L) 

35 CONTINUE 
30 CONTINUE 
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25 CONTINUE 
DO 36 1=1, K 

***IDIS = 1 FOR DISCRETE VARIABLES*** 

READ(10,*)IDIS(I) 

IF(IDIS(I).EQ.l) THEN 
♦♦♦DISCRETE FACTOR BY VARIABLE*** 

READ(10,*)RJ(D 
END IF 

36 CON TINUE 

♦♦♦INITIAL PENALTY FUNCTION ACCELERATION FACTOR*** 

ACCEL = 1.0 

IF(IOPT.EQ.l) THEN •. . mil* -iw . _ V- 

♦♦♦RANDOM SEARCH*** 

CALLRSEARCH(IDIS,ISEED,N,K,P,C,A,M,XKL,CI,AI,RJ,ACCEL,X) 

GO TO 1000 

ENDIF ' • . 

IF(IOPT.EQ.2)THEN 
♦♦♦HOOKE AND JEEVES*** 

CALL HJ(IDIS,N,K,P,C AM,XKL,CI,AI,RJ,ACCEL,X) 

ENDIF 

1000 CONTINUE 
10 CONTINUE 
STOP 
END 

SUBROUTINE RSEARCH(IDIS,ISEED,N,K,P,C,A,M,XKL,CI,AI,RJ,ACCEL,X) 
DIMENSION C(100),A(100, 100), X(100),M(100),CI(100,100) 

DIMENSION AI(100,100,100),XKL(100)^HIGH(100)^LOW(100) 

DIMENSION RJ(100),IDIS(100),XOPT(100),XEXP(100) 
DIMENSIONEPS1(100),EPS2(100),MULT(100) 

♦♦♦FRACTION OF INTERVAL REQUIRED AND CONFIDENCE LEVEL*** 
READ(10,*)FRS,XPRS 
♦♦♦NUMBER OF SEARCH POINTS*** 

WRrrE(6,*)’YOU WILL BE SEARCHING ’ .NPOINTS, ’POINTS’ 

DO 40 1=1, K 

♦♦♦LOWER AND UPPER BOUNDS BY VARIABLE*** 

READ( 1 0,*)XLOW(I),XHIGH(I) 

40 CONTINUE 

♦♦♦INITIALIZING NUMBER OF DISCRETE POINTS AND VARIABLES*** 

DPOINTS=1.0 
NDVAR=0 
DO 41 1=1, K 

IF(IDIS(I).EQ.l)TTffiN . ...... 

NDVAR=ND VAR+ 1 

♦♦♦CALCULATES TOTAL NUMBER OF FEASIBLE DISCRETE POINTS IN INTERVAL*** 
DPOINTS=(XHIGH(I)-XLOW(I))/RJ(I)*DPOINTS 
ENDIF 

41 CONTINUE 

***IF THE PROBLEM ISNT MIXED*** 
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105 IF(NDVAR.EQ.K)THEN 

♦♦♦IFTHE INTERVAL IS NOT DENSE IN DISCRETE FEASIBLE POINTS RELATIVE TO*** 
***THE NUMBER YOU WERE WILLING TO SEARCH ANYWAY, JUST SEARCH FEASI- 
BLE POINTS*** 

EF(DPOINTS.LE.NPOINTS)THEN 
DO 42 1=1 JC 

MULT(I)=IFDC(XLOW(I)/RJ(I))+ 1 
***LOWEST DISCRETE FEASIBLE POINT*** 

X(I)=MULT(I)*RJ(I) 

42 CONTINUE 
ICOUNT=0 
DO 44 J=1,K 

♦♦♦CONTINUE AS LONG AS DISCRETE POINTS ARE FEASIBLE*** 

47 IF(X(J).LE.XHIGH(J))THEN 

CALL OBJ (IDIS ,C,A,X,CI, AI,XKL,N,K,P,M,RJ,FUNC,ACCEL,XPF) 
IF(ICOUNT.EQ.O)THEN 
♦♦♦INITIALIZE OPTIMAL VALUES*** 

FUNCOPT=FUNC 

XPFOPT=XPF 


DO 43 1=1, K 
XOPT(I)=X(I) 

43 CONTINUE 
ENDIF 

ICOUNT=ICOUNT+l 
IF(FUNC.LT.FUNCOPT)THEN 
♦♦♦UPDATE OPTIMAL VALUES*** 

FUN COPT =FUN C 
XPFOPT=XPF 
DO 46 L=1,K 
XOPT(L)=X(L) 

46 CONTINUE 

ENDIF 

♦♦♦INCREMENT DISCRETE SEARCH POINTS*** 

MULT (J)=MULT ( J)+ 1 
X(J)=MULT(J)*RJ(J) 

GO TO 47 
ENDIF 

♦♦♦UPDATE OPTIMAL VALUES*** 

DO 48 1=1, K 
X(I)=XOPT(I) 

48 CONTINUE 

44 CONTINUE 

*** WRITE(6,*)FUNCOPT ,XPFOPT,(XOPT(I),I= 1 ,K) 

GO TO 99 

ENDIF 

ENDIF 

♦**IF THE PROBLEM IS MIXED OR CONTINUOUS OR FULLY DISCRETE WITH A DENSE 
♦♦♦COVERING OF FEASIBLE POINTS IN THE INTERVAL, PROCEED WITH STANDARD 


RANDOM 



An Employee-Owned Company 


148 


♦♦♦SEARCH*** 

D045 1=1 ^POINTS - • ’ 

DO 50 J=1,K 

♦♦♦CALCULATE RANDOM SEARCH POINT*** 

X(J)=XLOW(J)+RAN(ISEED)*(XHIGH(J)-XLOW(J)) 

50 CONTINUE 

CALLOBJ(IDIS,C,A,X,CI,AI,XKL,N,K,P,M,RJ,FUNC,ACCEL,XPF) 

IF(I.EQ.1)THEN 

♦♦♦INITIALIZE OPTIMAL VALUES*** 

FUNCOPT=FUNC 
XPFOPT=XPF 
DO 85 L=1,K 
XOPT(L)=X(L) 

85 CONTINUE 

ENDIF '-t • 

IF(FUNC.LT.FUNCOPT)THEN 
♦♦♦UPDATE OPTIMAL VALUES*** 

FUN COPT =FUN C 
XPFOPT=XPF 
DO 90 L=1,K 
XOPT(L)=X(L) 

90 CONTINUE 
ENDIF 

**♦ WRHE(6 *)XPFOPT 

♦♦♦DOES PENALTY CONVERGE?*** 

99 IF(XPFOPT.LE.O.OO 1 )THEN 
WRITE(11,*)’MIN. OBJ. FUNC. VALUE =’,FUNCOPT 

DO 95 L=1,K .r . 

WRITE(1 1 ,*)’X’,L, , =’,XOPT(L) 

95 CONTINUE 
GOTO 100 
ENDIF 

♦♦♦UPDATE PENALTY FUNCTION ACCELERATING FACTOR IF PENALTY DOESNT 
♦♦♦CONVERGE*** 

ACCEL=ACCEL* 10.0 
GOTO 105 

100 CONTINUE 
RETURN 
END 

SUBROUTINE HJ(IDIS,N,K,P,C,A,M,XKL,CI,AI,RJ,ACCEL,X) 

DIMENSION C(100),A(100, 100), X(100),M(100),CI(100,100) 

DIMENSION AI(100,100,100),XKL(100),XHIGH(100),XLOW(100) 

DIMENSION RJ(100),IDIS(100),XOPT(100),XEXP(100) 

DIMENSION EPS 1(100), EPS2( 100) 

IDIFF=0 
DO 110 1=1, K 

♦♦♦READ INITIAL POINT, INITIAL EXPLORATORY VALUES, AND FINAL EXPLOR- 


ATORY 
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READ(10 ,*)X(I),EPS 1(I)3PS2(D 
XOPT(D=X(I) 

110 CONTINUE 

CALLOBJ(IDIS,C,A,X,CI,AI,XKL,N,K,P,M,RJ,FUNC,ACCEL,XPF) 
*** WRITE(11,*)T,FUNC 
♦♦♦INITIALIZE OPTIMAL VALUES*** 

FUN COPT =FUN C 
XPFOPT=XPF 

*** WRITE( 1 1 ,*)FUNCOPT ,XPFOPT 
DO 115 1=1, K 

♦♦♦PERFORM EXPLORATORY SEARCH FROM BASE POINT*** 
XEXP(I)=X(I) 

x(i)=x(i)+EPsia) 

CALLOBJ(IDIS,C,A,X,CI,AI,XKL,N,K,P,M,RJ,FUNC,ACCEL,XPF) 
*** WRITE(1 1,*)’2\FUNC 

IF(FUNC.GT.FUNCOPT)THEN 
***GO IN OTHER DIRECTION*** 

X(I)=X(I)-2.0*EPS 1(1) 

*** DO 1134 KLM=1,K 

*** WRITE( 11,*)’ KLM= ’ ,KLM,X(KLM) 

*** 1134 CONTINUE 

CALLOBJ(roiS,CAX,CI,AI,XKL,N,K,P,M,RJ,FUNC,ACCEL,XPF) 
*** WRITE(1 1,*)’3\FUNC 

IF(FUNC.LT.FUNCOPT)THEN 
♦♦♦UPDATE OPTIMAL VALUES*** 

FUNCOPT=FUNC 

XOPT(I)=X(I) 

XPFOPT=XPF 

IDIFF=1 

*♦* DO 1125 JIK=1,K 

*** WRITE(1 1,*)’ 1 1 1 ’,FUNCOPTXOPT(JIK),XPFOPT 

*** 1125 CONTINUE 
GOTO 111 
END IF 
GOTO 111 
ENDIF 

♦♦♦UPDATE OPTIMAL VALUES*** 

XOPT(I)=X(I) 

XPFOPT=XPF 
FUN COPT =FUNC 
IDIFF=1 

*** DO 1126 JIK=1,K 

*** WRITE(11 *)’115’J 7 UNCOPT^OPT(JIK),XPFOPT 

*** 1126 CONTINUE 

111 CONTINUE 
115 CONTINUE 
135 IF(IDIFF.EQ. 1 )THEN 

DO 120 1=1, K 
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***DF NEW POINT IN EXPLORATION IS DIFFERENT FROM BASE POINT, PERFORM 
♦♦♦PATTERN SEARCH*** 

XEXP(I)=XEXP(I)+2.0*(XOPT(I)-XEXP(I)) 

X(I)=XEXP(I) 

120 CONTINUE 

ENDIF •• '• •' --- 

♦♦♦OTHERWISE, REDUCE EXPLORATORY VALUES (EPSILONS)*** 
IF(IDIFF.EQ.O)THEN 
GOTO 140 
ENDIF 
ID1FF=0 

CALLOBJ(IDIS,C,A,X,CI,AI,XKL,N,K,P,M,RJ,FUNC,ACCEL,XPF) 

*** WRITE(1 1,*)’4\FUNC " 

IF(FUNC.LT.FUNCOPT)THEN 
♦♦♦UPDATE OPTIMAL VALUES*** 

FUNCOPT=FUNC - ^ - • 

XPFOPT=XPF 
DO 121 I=1,K 

XOPT(I)=X(I) - ----- - 

*** DO 1127 JIK=1,K 

*** WRITE( 1 1 ,*) T 2 1 ’ JFUNCOPT ,XOPT (JIK),XPFOPT 

*** 1127 CONTINUE 

121 CONTINUE 
175 DO 125 1=1, K 

♦♦♦PERFORM (NEW) EXPLORATION*** 

X(I)=X(I)+EPS1(I) " 

CALLOBJ(TOIS,C,A,X,CI,AI,XKL,N,K,P,M,RJ,FUNC,ACCEL,XPF) 

*** WRITE(1 1,*)’5’,FUNC 

IF(FUNC.GT.FUNCOPT)THEN 
***GO IN OTHER DIRECTION*** 

X(I)=X(I)-2.0*EPS1(I) 

*** DO 1136 KLM=1,K 

*** WRITE(1 1,*)’KLM=’,KLM,X(KLM) 

*** 1136 CONTINUE 

CALL OBJ(roiS,CA,X,CI,AI,XKL,N,K,P,M,RJ,FUNC,ACCEL,XPF) 

*** WRITE(ll,*)’6’JFyNC 

IF(FUNC.LE.FUNCOPT)THEN 
♦♦♦UPDATE OPTIMAL VALUES*** 


FUN COPT -FUNC 
XOPT(I)=Xa) 

XPFOPT=XPF 
*** DO 1 128 JIK=1,K 

*** WRITE(ll,*)T75 , EUNCOPT^OPT(JIK) v XPFOPT 

*** 1128 CONTINUE 


IDIFF=1 
ENDIF 
GO TO 125 
ENDIF 

♦♦♦UPDATE OPTIMAL VALUES*** 
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XOPT(I)=X(I) 

IDIFF=1 

FUN COPT =FUN C 
XPFOPT=XPF 
*** DO 1129 JIK=1,K 

*** WRITE( 1 1 , *) ’ 1 40’ ,FUNCOPT,XOPT(JIK),XPFOPT 

*** 1129 CONTINUE 
125 CONTINUE 
140 IF(IDIFF.EQ.O)THEN 
DO 130 1=1, K 

♦♦♦REDUCE EXPLORATORY VALUES WHEN NO IMPROVEMENT IS MADE IN 
EXPLORATION*** 

EPS 1 (I)=EPS 1 (I)/2.0 
IF(EPS 1 (I).LT.EPS2(I))THEN 

♦♦♦CHECK ENDING CONDITION BASED ON EXPLORATORY VALUES*** 

GOTO 150 
END IF 

130 CONTINUE 

***GO EXPLORE SOME MORE*** 

GO TO 175 
END IF 
END IF 

DO 137 KIM=1,K 

♦♦♦RETAIN PREVIOUS BASE POINTS FOR FUTURE PATTERN MOVES*** 
X(KIM)=XOPT(KIM) 

137 CONTINUE 

IF(IDIFF,EQ.0)THEN 
GO TO 140 
END IF 

***MAKE PATTERN MOVE*** 

GOTO 135 
150 CONTINUE 

♦♦♦CHECK PENALTY VALUE FOR CONVERGENCE*** 

IF(XPFOPT.LE.0.001)THEN 
DO 1130 JIK=1,K 

WRITE(1 1,*)’MIN. OBJ. FUNC. VALUE =\FUNCOPT 
WRITE(11,*)’PENALTY = ’.XPFOPT 
WRITE(1 1,*)’ ACCELERATION FACTOR = ’ACCEL 
1130 CONTINUE 
DO 170 L=1,K 

WRITEG 1 *)’X’,L,’=’ AOPT(L) 

170 CONTINUE 
GO TO 200 
ENDIF 

♦♦♦UPDATE PENALTY FUNCTION ACCELERATION FACTOR IF CONVERGENCE IS NOT 
***ACHIEVED*** 

ACCEL= ACCEL* 10.0 
GOTO 110 
200 CONTINUE 
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RETURN 
END 

SUBROUTINE OBJ(IDIS,CA^,CI,AI,XKL^,K^M,RJ,FUNC,ACCEL,XPF) 
DIMENSION C(100)A(100,100),X(100),M(100),a(100,100) 

DIMENSION AI(100,100,100)^KL(100)^ffiGH(100)^LOW(100) 
DIMENSION RJ(100),IDIS(100),XOPT(100),XEXP(100) 

DIMENSION EPS 1 ( 1 00),EPS2( 1 00) 

FUNC=0.0 
DO 55 U=1,N 

♦♦♦INITIALIZE PRODUCT VALUES*** 

XPROD=1.0 

DO60UK=l,K r _ - - 

****** WRITE(6,*)X(UK)A(U,UK) 

***ZERO EXPONENTS GIVE PRODUCT VALUES OF 1.0*** 

IF(AB S ( A(U 4JK)).LE.0.000000 1 )THEN 
GO TO 60 ' 

ENDIF 

***ZERO VARIABLE VALUES GIVE PRODUCT VALUES OF ZERO*** , „ 
IF(ABS(X(UK)).LE.0.0000001)THEN 
XPROD=0.0 
GO TO 60 
ENDIF 

♦♦♦COMPUTERS DONT RAISE NEG. VALUES TO EXPONENTS*** 
IF(X(UK).LT.0.0)THEN 
X(UK)=AB S (X(I JK)) 

ENDIF 

♦♦♦COMPUTE PRODUCT VALUES FOR OBJ. FUNCTION*** 
XPROD=XPROD*X(UK) * * A(U,I JK) 

60 CONTINUE 

♦♦♦COMPUTE ORIGINAL OBJECTIVE FUNCTION VALUE*** 
FUNC=FUNC+C(U)*XPROD 
55 CONTINUE 
XPF=0.0 
DO 65 U=1,P 

♦♦♦INITIALIZE CONSTRAINT SUMS*** 

CONSUM=0.0 
DO 70 LJ=1,M(U) 

♦♦♦INITIALIZE CONSTRAINT PRODUCTS*** 

CONPROD=1.0 
DO 75 MJ=1,K 

♦♦♦UPDATE CONSTRAINT PRODUCTS*** 

CONPROD=CONPROD*X(MJ)**AI(U,MJ,U) 

75 CONTINUE 
***UPDATE CONSTRAINT SUMS*** 

CON SUM=CON SUM+CI(U,U) * CONPROD 
70 CONTINUE 

♦♦♦COMPUTE PENALTY AND FUNCTION FOR <= CONSTRAINTS*** 
IF((CONSU M-XKL(U )).GT.0.0)THEN 
XPF=XPF+ACCEL* (CONSUM-XKL(IJ» * *2.0 
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FUN C=FUN C+XPF 
ENDIF 

65 CONTINUE 
DO 80 U=1,K 

♦♦♦COMPUTE PENALTY AND FUNCTION FOR DISCRETE CONSTRAINTS*** 
IF(IDIS(U).EQ.1)THEN 

XPF1 =ACCEL*(ABS (X(U)/RJ(U)*IFIX(X(U)/RJ (U)+0.5)))**0.5 

XPF=XPF+XPF1 

FUNC=FUN C+XPF 1 

ENDIF 

80 CONTINUE 
RETURN 
END 
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APPENDIX B. IMPACT 10 SOURCE CODE LISTING 

DIMENS ION XPV(IOO), SOLAR(l 188), XPSIV(125), XMETIV(IOO) 

DIMENSION XDEBOLDIV ( 1 00) 

DIMENSION ISWITCH(IO) 

character BUMPER_NAME(50)*14,WALL_NAME(50)*1 1 

CHARACTER BUMPER_MAT_NAME*40 , WALL_MAT_NAME*40 

CHARACTER BUMPER JTYPE_NAME*40,WALL_TYPE_NAME*40 

CHARACTER SHAPE*40 

character Outdir*40 

character line 1*80 

character line2*80 

data cdate / ’Run_Date 7 

data crime / ’RunJTime 7 

call gettim (ihr.imin, isec, ilOOth) 

call getdat (iyr, imon, iday) 

OPEN(UNIT=27,STATUS= , old , ^CCESS=’SEQUENTIAL^ 

+ FILE=’confIg.pgm’) 
read( 27 y 231 2)outdir 
read(27,23 12)outdir 
2312 format(A40) 
close(27) 

OPEN(UNIT=23,STATUS=’OLD\ACCESS=’SEQUENTIAL\FILE=’CRAFT.INP’) 
OPEN(UNIT=26,STATUS=’OLD’,ACCESS=’SEQUENTIAL’,FILE=’GEOMETRY.INP’) 
UK - INDEX(OUTDIR, ’ ’) -1 

OPEPr(UNIT=27,STATUS= , unknown’ACCESS= , SEQUENTIAL’, 

+ HLE= outdir(l:UK) // ’Z9AAAAAJ.PGM’) 

OPEN (UNIT=2 8 ,STATU S = ’ unknown ’ ,ACCES S = ’ SEQUENTIAL’ , 

+ HLE=’PROJECT.OUT) 

OPEN(UNIT=29,STATUS=’unknown’ACCESS=’SEQUENTIAL’, 

+ FILE=’results.dat’) 

open(unit=33,status= ’old ’ ,access= ’ sequential ’ , 

+ file= ’project, hdr’) 

write(28, ’(lx, A10,lx, 12.2, 1H:, 12.2, 1H:, 12.2, 1H., 

+12.2)’) crime, ihr, imin, isec, ilOOth 

write(28, ’(lx, A10,lx, 12.2, 1H-, 12.2, 1H-, 14.2)’) 

+ cdate, imon, iday, iyr 
do 6008 i = 1,6 
read(33 ,6007)line 1 
write(28,6007)line 1 
6008 continue 
close(unit=33) 

write(*,*)’ IMPACT10V - SAIC / Huntsville’ 

write(*,*)’ ’ 

write(*,*) 

write(*,*) ’Status - Initializing Files’ 
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C READ Static Data Files 
C 

OPEN(UNIT=14,STATUS=’OLD’,ACCESS=’SEQUENTIAL\- 

FILE=’FLUXFAC.DAT’) 

DO 222 KI=1,101 
READ(14,*)JI,XPSIV(JI) 

222 CONTINUE 
CLOSE (UNIT = 14) 

OPEN(UNIT=14,STATUS=’OLD’,ACCESS=’SEQUENTIAL\FILE=’SOLARl.FLX’) 
DO 223 KI=1,1188 
READ(14,*)SOLAR(KI) 

223 CONTINUE 
CLOSE (UNIT = 14) 

OPEN(UNIT=14,STATUS=’OLD’,ACCESS=’SEQUENTIAL’,FILE=’METVEL.INP’) 
DO 224 KI= 1,72 
RE AD( 1 4,*)XMETTV (KI) 

224 CONTINUE 
CLOSE (UNIT = 14) 

OPEN(UNIT=14,STATUS=’OLD\ACCESS=’SEQUENTIAL\ 

+ HLE=’DEBOLDVE.DAT’) 

DO 225 KI=1,16 

READ(14,*)IV,XDEBOLDIV(IV) 

225 CONTINUE 
CLOSE (UNIT = 14) 

c 

c Read PSDOC controlling switches and set impactlO variables accordingly 
c 

OPEN(UNIT=14,STATUS= , OLD\ACCESS=’SEQUENTIAL\ 

+ FILE=’SWITCH.INP’) 

DO 226 KI=1,10 
READ(14,*)ISWITCH(KI) 

226 CONTINUE 
CLOSE (UNIT = 14) 

NENVIRON = ISWITCH(l) 

IBUMPER_TYPE=ISWITCH(2) 

IBUMPER_MATERIAL=ISWITCH(3) 

IWALL.TYPE = IS WITCH (4) 

IWALL_MATERIAL=ISWITCH(5) 

isw6 = iswitch(6) 
isw7 = iswitch( 7) 
isw8 = iswitch(8) 

C NOTE: CURRENTLY WE ARE NOT MAKING USE OF ISWITCH(9) 
IGRAPH_TYPE=IS WITCH ( 1 0) 


C Open appropriate files for environment depending on switch 
c settings. 
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IF (NENVIRON.EQ. 1 ) THEN 

OPEN(UNIT=20,STATUS=’OLD’,ACCESS=’SEQUENTIAL\ 
+ FILE=’NEWDEBRUNP’) 

END IF 

IF (NENVIRON.EQ.2) THEN 

OPEN (UNIT=20, STATUS = ’ OLD ’ , ACCES S = ’ SEQUENTIAL’, 
+ FILE=’OLDDEBRI.INP’) 

END IF 

IF (NENVIRON.EQ.3) THEN 

OPEN(UNIT=2(),STATUS=’OLD\ACCESS=’SEQUENTIAL\ 
+ FILE=’NE_MET.INP’) 

END IF 

C Open appropriate files and read data from bumper database if 
c the table data is used rather than the single material (parametric) 
settings. __ ■ : 

C- 

BUMPER DATABASE 


IF (IBUMPER_TYPE.EQ. 1) THEN 
BUMPER_TYPE_NAME=’Bumper Material Database’ 

IF (IBUMPER_MATERIAL.EQ. 1 ) THEN 
BUMPER_MAT_NAME=’ Aluminum Alloy’ 
OPEN(UNIT=30,STATUS=’OLD’,ACCESS=’SEQUENTIAL\ 

+ FILE=’ALBUMP.INP’) 

OPEN(UNIT=32,STATUS=’OLD’,ACCESS=’SEQUENTIAL\ 

+ FILE= , ALBM.TBL’) 

ELSE IF (EBUMPER_material.EQ.2)THEN 
BUMPER_MAT_NAME=’Titanium Alloy’ v . . 

OPEN(UNIT=30,STATUS=’OLD’ ) ACCESS= 1 SEQU£NTIAL’, 

+ FILE=’TIBUMP.INP’) 

OPEN(UNIT=32,STATUS=’OLD ’ , ACCES S=’SEQUENTIAL’, 

+ FRE=”nTBM.TBL’) 

ELSE IF (ffiUMPER_MATERIAL.EQ.3)THEN 
BUMPER_MAT_NAME=’Steel Alloy’ 
OPEN(UNIT=30,STATUS=’OLD\ACCESS=’SEQUENTCAL\ 

+ FILE=’STBUMP.INP’) 

OPEN(UNIT=32,STATUS=’OLD’,ACCESS=’SEQUENTIAL’, 

+ FILE=’STBM.TBL’) ::i: ^ 

ELSE IF (IBUMPER.MATERIAL.EQ.4) THEN 
BUMPER_MAT_NAME=’Inconel Alloy’ 
OPEN(UNIT=30,STATUS=’OLD’,ACCESS=’SEQUENTIAL’, 

+ FILE=TNBUMP.INP’) • 

OPEN (UNTT =32,STATU S =’ OLD ’ , ACCES S= ’ SEQUENTIAL’ , 

+ FTLE=’INCBM.TBL’) 

ELSE IF (IB UMPER_MATERIAL. EQ.5 )THEN 
BUMPER_MAT_NAME=’Graphite Alloy’ 
OPEN(UNIT=30,STATUS=’OLD\ACCESS=’SEQUENTIAL\ 

+ FILE=’GRBUMP.INP’) ^ * A |g 
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OPEN(UNIT=32,STATUS=’OLD\ACCESS=’SEQUENTIAL\ 
+ FILE=’GRALBM.TBL’) 

ELSE 

write(*,*)’Input Error from PSDOC interface - ’ 
write(*,*) ’Program Terminated with internal error.’ 
write(*,*) ’Bad Ibumper_material switch, 
goto 1 1 
END IF 

C Read Number of bumpers in selected table 
Read(30,31)nbump 
31 format(IlO) 

c Read Table files for Material names 
c Skip 7 line header first 

DO 5010 KI= 1,7 
READ(32,*) 

5010 CONTINUE 

DO 627 KI=l,nbump 

READ(32,6000,END=628)BUMPER_NAME(KI) 

627 CONTINUE 
6000 format(al4) 

628 CLOSE (UNIT = 32) 


c For parametric settings, open bumper.inp using same file handle as table 
c setting. This will allow us to use the same code regardless of the method 
c chosen. 

C 

C SINGLE BUMPER MATERIAL 

C 

else IF (IBUMPER_TYPE.EQ.2) THEN 

BUMPER_TYPE_NAME=’ Single Bumper Material’ 
OPEN(UNIT=30,STATUS=’OLD\ACCESS=’SEQUENTIAL\ 

+ FILE=’BUMPER.INP’) 

NBUMP = 1 
else 

write(*,*) ’Input Error from PSDOC interface - ’ 
write(*,*)’Program Terminated with internal error.’ 
write(*,*)’Bad Ibumper_type switch, 
goto 1 1 
END IF 

C Open appropriate files and read data from wall database if 
c the table data is used rather than the single material (parametric) 
c settings. 

C 

£*************** 

C SELECT WALL MATERIAL (SINGLE /DATABASE) 

************** 
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C 

C 

IF (IWALL.TYPE.EQ. 1 )THEN 
WALL_TYPE_NAME=’Wall Material Database’ 

IF (IWALL_MATERIAL.EQ. 1 )THEN 
WALL_MAT NAME=’Aluminum Alloys’ 
OPEN(UNIT=35,STATUS=’OLD\ACCESS=’SEQUENTIAL\ 
+ FILE=’ALWALL.INP’) 

OPEN(UNIT=39,STATUS=’OLD\ACCESS=’SEQUENTIAL\ 
+ FTLE=’alwall.tbr) 

else IF (IWALL_MATERIAL.EQ.2)THEN 
WALL_MAT_NAME=’Advanced Launch System’ 
OPEN(UNIT=35,STATUS=’OLD’,ACCESS=’SEQUENTIAL’, 
+ FILE=’ALSWALL.INP’) 

OPEN (UNIT =39, STATU S = ’OLD ’ ,ACCES S=’ SEQUENTIAL \ 
+ FILE=’ALSWALL.TBL’) 

else 

write(*,*) ’Input Error from PSDOC interface - ’ 
write(*,*) ’Program Terminated with internal erro r. ’ 
write(*,*)’Bad Iwall_material switch. ’ 
goto 1 1 
end if 

C Read Number of walls in selected table 
Read(35,31)nwall 

DO 5015 KI=1,7 

READ(39,*) ‘ 

5015 CONTINUE 

DO 637 KI=1,NWALL 

READ(39,6002,END=638)WALL_NAME(KI) 

637 CONTINUE 

638 CLOSE (UNIT = 39) 

6002 format(all) 


C SINGLE WALL MATERIAL 

else IF (IWALL TYPE.EQ.2)THEN 
WALL_TYPE_NAME=’Single Wall Material’ 
OPEN(UNIT=35,STATUS=’OLD’,ACCESS=’SEQUENTIAL’, 
+ FILE=’WALL.INP’) 

NWALL = 1 
else 

write(*,*)Tnput Error from PSDOC interface - ’ 
write(*,*)’Program Terminated with internal error.’ 
write(*,*)’Bad Iwall_type switch. ’ 
goto 1 1 
END IF 
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c if using Table option, open new table 1. out file and write headers 
if ((ibumper_type.eq.l).or.(Iwall_type.eq.l)) then 
OPEN(UNIT=37,STATUS=’unknown’,ACCESS=’SEQlJENTIAL’, 
+ FILE=TABLE1 .OUT’) 

C Read headers into table 1. out file 

open(unit=33,status=’old’,access=’sequential’, 

+ file=’tablel.hdr’) 

write(37, ’(lx, A 10, lx, 12.2, 1H:, 12.2, 1H:, 12.2, 1H., 

+12.2)’) ctime, ihr, imin, isec, ilOOth 

write(37, ’(lx, A10,lx, 12.2, 1H-, 12.2, 1H-, 14.2)’) 

+ cdate, imon, iday, iyr 
do 6005 i = 1,4 
read(33,6007)line 1 
write(37,6007)line 1 

6005 continue 
6007 format(a80) 

read(33,6007)linel 

line2 = linel(l:12) // bumper_mat_name // linel(53:80) 

write(37,6007)line2 

read(33,6007)line 1 

line2 = linel(l:12) // wall_mat_name // linel(53:80) 
write(37,6007)line2 
do 6006 i=l, 7 
read(33,6007)line 1 
write(37,6007)line 1 

6006 continue 
close(unit=33) 

end if 
C 

i************* 

C READ GEOMETRICAL SHAPE 

£»* a|t * * * * * * * * * * * * * 

C 

IF (ISW6.EQ.1) THEN 
SHAPE = ’Cylinder’ 
else 

write(*,*)’Input Error from PSDOC interface - ’ 
write(*,*) ’Program Terminated with internal error.’ 
write(*,*)’Bad Isw6 switch, 
goto 1 1 
END IF 


* * * * * * *** * * * * 

C READ IMPACT MODEL 

Q*4^************ 

c 
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IF (ISW7.EQ.1) THEN 
IMPACT_MODEL=’Single Impact Mo del’ 
else IF (ISW7.EQ.2) THEN 
IMPACT_MODEL= ’Three Impact Regions’ 
else 

write(*,*) ’Input Error from PSDOC interface - ’ 
write(*,*)’Program Terminated with internal error.’ 
write(*,*)’Bad Isw7 switch. ’ 
goto 1 1 
end if 


IF (ISW8.EQ.1) THEN 
DEBRIS_ENVIRONMENT=’Boeing Model’ 

NCODE = 2 
else 

write(*,*)Tnput Error from PSDOC interface - ’ 
write(*,*) ’Program Terminated with internal error.’ 
write(*,*)’Bad Isw8 switch, 
goto 1 1 
END IF 
C- 

C SELECT GRAFTOOL OUTPUT FILE FORMAT 

e 

C Write the definitions and the headings in the temporary output 
C file ( z9aaaaaj.pgm ) in the c:\psdoc\mput directory, 
c This format is critical to PSDOC front-end...don’t change! 

IF (IGRAPH_TYPE.EQ. 1 )THEN 
WRTTE(27,*)7 The Output Variables Are Defined As: ’ 

WRITE(27 ,*)’/ ’ 

WRITE(27,*)7 Column 1: run = Run-#’ 

WRITE(27,*)7 Column 2: T1 = Optimal Bumper Thickness’ 
WRITE(27,*)7 Column 3: T2 = Optimal Wall Thickness’ 
WRITE(27,*)7 Column 4: OBMPUA = Optimal Bumper Mass’// 
+ ’ Per Unit Area’ 


+ 

+ 


+ 


+ 


WRITE(27,*)7 Column 5: OWMPUA = Optimal Wall Mass’// 


’ Per Unit Area’ 
WRITE(27,*)7 Column 6: 
//’Unit Area’ 
WRITE(27,*)7 Column 7 : 
WRITE(27 ,*)’/ Column 8 : 
WRITE(27*)7 Column 9 : 
WRITE(27 *)’/ Column 10 : 
// ’ Diameter’ 

WRITE(27 *)’/ Column 11: 
WRITE(27,*)7 Column 12: 
// ’ Rate’ 


WT = Minimum System Mass Per’ 

WTCMC = Mnimum CMC Weight’ 
OBR = Optimal Bumper Ratio’ ~ 
OWR = Optimal Wall Ratio’ 

D = Critical Design Projectile’ 

RHOP = Projectile Density’ 
XGROWTH = Space Debris Growth’ 
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WRTTE(27,*)7 Column 13: IMONTH1 = Initial Operation Month’ 
WRITE(27 ,*)’/ Column 14: I YEAR 1 = Initial Operation Year’ 
WRITE(27,*)7 Column 15: IMONTH2 = Final Operation Month’ 
WRITE(27,*)7 Column 16: IYEAR2 = Final Operation Year’ 
WRITE(27 ,*)’/ Column 17: ALT = Spacecraft Orbital’// 

+ ’ Altitude’ 

WRITE(27,*)7 Column 18: XINCL = Spacecraft Orbital’// 

+ ’ Inclination’ 

WRTTE(27 ,*)’/ Column 19: XPO = Spacecraft Probability’// 

+ ’ Of No Penetration’ 


WRITE(27,*)7 Column 20: AREAK = Spacecraft Exposed’// 
+ ’ Area’ 


+ 


+ 

+ 


+ 

776 


C 

C 


+ 


+ 

+ 

+ 


+ 


WRITE(27,*)7 Column 21 : S = Spacecraft Bumper / Wall’// 

’ Separation’ 

WRITE(27 ,*)’/ ’ 

WRITE(27 ,*)’/ ’ 

WRITE(27 ,*)’/ ’ 

WRITE(27,*)7 Valid X- Columns Are: ’ 

WRITE(27,*)7 11,16,17,18,19,20 ’ 

WRITE(27,*)7Valid Y- Columns Are: ’ 

WRITE(27,*)7 2,3,4, 5, 6, 7,8, 9 ’ 

WRUE(27 ,*)’/ ’ 

WRITE(27,*)7 ’ 

WRITE(27,*)7 ’ 

WRITE(27,776)’/RUN-# ’,’ T1 ’,’ T2 ’,’ OBMPUA ’,’OWMPUA’, 
’WT’.’WTCMC’,’ OBR ’,’ OWR ’,’ D ’,’ RHOP ’,’XGROWTH’, 
’IMONTH1’,’ IYEAR1’,’ IMONTH2’,’IYEAR2’ 

,’ ALT ’,’ XINCL ’,’ XPO ’,’ AREAK ’,’ S ’ 

FORMAT(21(A12,lX)) 

C 

WRITE THE DEFINITIONS AND THE HEADINGS IN THE OUTPUT 
FILE ( RESULTS. OUT ) IN THE LOTUS (123) FORMAT 
else IF (IGRAPH_TYPE.EQ.2)THEN 
WRITE(27,*)’ " ’,’ The Output Variables Are’// 

’ Defined As: ’,’ " ’ 

WRITE(27 ,*) 

WRITE(27 ,*)’ " ’,’ T1 = Optimal Bumper Thickness’,’ " ’ 

WRITE(27 ,*)’ " ’,’ T2 = Optimal Wall Thickness’,’ " ’ 

WRITE(27 *)’ " ’,’ OBMPUA = Optimal Bumper Mass Per Unit ’ 

// ’Area’ ’ " ’ 

WRITE(27,*)’ " ’,’ OWMPUA = Optimal Wall Mass Per 7/ 

’Unit Area’,’ " ’ 

WRITE(27,*)’ " ’,’ WT = Minimum System Mass Per Unit 7/ 

’Area’ ’ " ’ 

WRITE(27,*)’ " ’,’ WTCMC = Minimum CMC Weight’,’ " ’ 
WRITE(27,*)’ " ’,’ OBR = Optimal Bumper RATIO’,’ " ’ 
WRITE(27,*)’ " ’,’ OWR = Optimal Wall RATIO’ ’ " ’ 

WRITE(27 *)’ " ’,’ D = Critical Design Projectile ’ 

’Diameter’,’ " ’ 
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WRTTE(27 ,*)’ " ’,’ RHOP = Projectile Density’,’ " ’ 

WRTTE(27,*)’ " ’,’ XGROWTH = Space Debris Growth Rate’,’ " ’ 
WRITE(27,*)’ " ’,’ IMONTH1 = Initial Operation Month’ 

WRITE(27,*)’ " ’,’ IYEAR1 = Initial Operation Year’ 

WRITE(27 ,*)’ " ’,’ IMONTH2 = Final Operation Month’ 

WRITE(27 *)’ M ’,’ IYEAR2 = Final Operation Year’ 

WRITE(27,*)’ " ’,’ ALT = Spacecraft Orbital Altitude’,’ " ’ 

WRTTE(27,*)’ H ’,’ XINCL = Spacecraft Orbital ’// 

+ ’Inclination’,’"’ ... ... 

WRITE(27,*)’ ’’ ’,’ XPO = Spacecraft Probability Of No ’// 

+ ’Penetration’,’ " ’ 

WRTTE(27 .*)’ " ’,’ AREAK = Spacecraft Exposed Area’,’ " ’ 

WRITE(27,*)’ " S = Spacecraft Bumper / Wall ’// 

+ ’Separation’,’ " ’ 

WRTTE(27,*)’ " " ’ 

WRTTE(27,*)’ " " ’ 

WRITE(27,*)’ " " ’ 

WRITE(27,*)’/ Valid X- Columns Are: ’ 

WRITE(27*)7 1 1,16,17, 18,19, 20 ’ 

WRITE(27*)7 Valid Y- Columns Are: ’ 

WRTTE(27,*)7 2,3,4,5,6,7,8,9 ’ 

WRITE(27,*) 

WRITE(27,*) 

WRTTE(27,*)’ 

WRITE(27,777)’ "RUN-#" ’,’"T1 " ’,’ "T2"’,’ "OBMPUA"’,’"OWMPUA , ”, 

+ ”*WT"’, ’ "WTCMC" ’ , ’ "OBR" ’ , ”’OWR" ’ , ’ "D" ’, ’ "RHOP" ’, ’"XGROWTH" ’ , 
+ ’ "IMONTH 1 ,n , ’ "IYEAR 1" ’ , ’ "IMONTH2" ’ , ’ "IYE AR2" ’, ’ "ALT" ’ , 

+ ’"XINCL"’, ’"XPO"’, ’"AREAK"’,’ ,, S’” 

777 FORMAT(2 1 (A 1 2, 1 X)) 

else 

write(*,*)’Input Error from PSDOC interface - ’ . 
write(*,*)’Program Terminated with internal error.’ 
write(*,*)’Bad Igraph switch, 
goto 1 1 
END IF 

C CALCULTE THE NUMBER OF MATERIALS " NMATS " 

NMATS = NBUMP * NWALL 


C IN HERE WE ARE INITIALIZING THE COUNTER VARIABLE "1" 


1=1 
iiii = 0 
write(*,*) 
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20 CONTINUE 

mi = nil 4 - 1 
write(*,21)iiii 

21 fonnat(’+Status - Running Case: ’,14) 
c 

c Reset files 20,23,26 (newdebri.inp, olddebri.inp, 
c ne_met.inp, craft.inp, and geometry.inp) 
c when running in table mode. If parametric, always rewind 30, 35 
c (bumper.inp, and wall.inp) for proper execution, 
if (ibumper_type.eq.l) then 
rewind(20) 
rewind(23) 
rewind(26) 

else if (ibumper_type.eq.2) then 
rewind(30) 
end if 

if (iwall_type.eq.2) then 
rewind(35) 
end if 

SEEDVAL = 73 
call Seed (SeedVal) 
tl =0 
t2 * 0 
obmpua = 0 
owmpua = 0 
wt =0 
wtcmc =0 
obr =0 

C 

IF (NCODE.EQ.1) then 
C NYSMITH 

C PROJECTILE DIAMETER IN CM **** READ(10,*) 
READ(22,*,end=l 1)D 

C BUMPER /Wall SEPARATION **** READ(10,*)H 
READ(23,*,end=l 1)S 
C **** READ(10,*)RHOl’ 

READ (24,*,end=l 1)RH01 
C **** RE AD( 1 0,*)RHO2 ’ 

READ (25,*,end=l 1)RH02 
C **** RE AD( 1 0,*)CMCLEN 

READ (26,*,end=l 1)CMCLEN 
C **** READ( 1 0,*)CMCRAD 

READ (26,*,end=l 1)CMCRAD 

c WRITE(11,*)’ NYSMITH’ 
c WRITE(11,*) 
c WRITE(11,*)’ INPUT’ 
c WRITE(11,*) 
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c WRITE(11 *)’ Projectile Diameter In CM = ’ 

c WRITE(1 1,*)’ Bumper/Wall Separation In CM = ’ ,H 

c WRITE(11 ,*)’ Bumper/Wall Separation In CM = \S 

c WRITE(11,*) 

TIT = 0 
T2T = 0 


DO 26 j=l,16 Z ^ M 
V = FLOAT(J) 

C CALL NYSMITH(VAH,RHO 1 ,RH02,T1 ,T2,WT,WTCMC) 

CALL NYSMITH(V, D, H, RHOl, RH02, Tl, T2) 

TIT ■ TIT + Tl * XPV(J) 

T2T = T2T + T2 * XPV(J) 

26 CONTINUE 


c 

c 

c 

c 

c 

c 

c 

c 


c 


T1-T1T 
T2 = T2T 

WT = RHOl * Tl + RH02 * T2 
R12 « CMCRAD 
R22 = R12 + T2 
Rll =R22 + H 

R21 =R11 +T1 — : 

WTCMC = 3.1416 * (CMCLEN / 1000.0) 

WTCMC=WTCMC*(RH01*(R21**2.-R11**2)+RH02*(R22**2.-R12**2.)) 
c WRITEd 1,*)’ OUTPUT’ 

WRITE(11,*) 

WRITE(1 1,*)’ Bumper Thickness = \T1,’CM’ 

WRITE(1 1 ,*)’ Wall Thickness = ’,T2,’CM’ 

WRITEd 1,*)’ Minimum Weight = ’/WT.’GM/Square CM’ 

WRITE(1 1 *)’ CMC Minimum Weight = \WTCMC,’KG’ 

WRITEd 1,*) 

WRITEd 1,*) 

WRITEd 1,*) 

C 

else IF (NCODE.EQ.2) then 
BOEING 

C NENVIRON = 1 => EARTH ORBITAL SPACE DEBRIS (NEW) 
IF(NENVIRON .EQ. 1) THEN 
READ(20 *,end=l l)XGROWTH 
READ(20,*,end=l l)x 
IMONTH1 = x 

READ(20,*,end=l l)x 
IYEAR1 =x 

READ(20,*,end=l l)x 
IMONTH2 = x 

READ(20,*,end=l l)x 
IYEAR2 = x 
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READ(20,*,end= 1 1 )ALT 
READ(20,*,end=l 1)XINCL 
READ(20,*,end=l 1)XP0 
READ(23,*,end=l 1)AREAK 

INCL = IFIX(XTNCL + .5) 

XPSI ■ XPSIV(INCL) 

CALL DEBRIS(XGR0WTH,S0LAR,XPSI,IM0NTH1 JYEARl JM0NTH2, 
+ IYEAR2ALT,XINCL^CP0,AREAK > D,XPV,IVMAX) 

RHOP = 2.8 
IF(D .GT. 1.0) THEN 
RHOP = 2.8/(D**0.74) 

END IF 

else if (NENVIRON .EQ. 2) THEN 

C NENVIRON = 2 ==> EARTH ORBITAL SPACE DEBRIS (OLD) 


READ(20,*,end=ll)T 
READ(20,*,end=l 1)XP0 
READ(23,*,end=l 1)AREAK 

CALL DEBRISOLD(T, XPO, AREAK, D) 
RHOP = 2.8 
DO 555 KIJK=1,16 
XPV(KUK) = XDEBOLDIV (KIJK) 

555 CONTINUE 

else if (nenviron.eq.3) then 
C NENVIRON = 3 — > Near Earth Meteoroid 

READ(23,*,end=l 1) AREAK 
READ(20,*,end=l 1)T 
READ(20 ,*,end=l 1)ALT 
READ(20*,end=ll)XP0 


DENS = .5 • 

CALL METEOROID(AREAK, T, XPO, ALT, DENS, D, L) 

RHOP = DENS 
IVMAX = 72 


DO 544 KUK=1,72 
XPV(KUK) = XMETTV (KIJK) 
544 CONTINUE 


end if 

READ(23,*,end=l 1)S 
READ(26,*,end= 1 1)CMCLEN 
READ(26,*,end=l 1)CMCRAD 
READ(30,*,end=l 1)RH01 
c RHOl = C_RH01(K) 
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READ(30,*,end=l 1)SY1 
c SY1 = C_SY1(K) 

READ(30,*,end=l 1)E1 
c El = C_E1(K) 

READ(35 ,*,end=l 1)RH02 
c RH02 = C_RH02(K) 

READ(35,*,end=l 1)XL2 
c XL2 = C_XL2(K) 

READ(35,*,end=l 1)SY2 
c SY2 = C_SY2(K) 


XN = .85 

SY1 = SY1 * 144000.0 
SY2 = SY2 * 144000.0 
El* El * 6.880285E+10 
TIT = 0.0 
T2T = 0.0 

DO 36 J=1,IVMAX 
V = FLOAT(J) 

EF(XINCL .GT. 40.0)THEN 
THETA = ACOS(-1.0 * V /IVMAX) - 1.57 
IF(THETA .GT. 1.57)THEN 
THETA = 1.57 
END IF 
ELSE 

THETA = ACOS(-1.0 * V / 15.4) - 1.57 
END IF 

C 3676 CALL BOEING(V .D.RHOPRHO 1 ,RH02,S,XL2,S Y 1 ,S Y2.THETA, 
C + XN.E1 .CMCRAD.Tl ,T2,WT,WTCMC) 


3676 CALL B0EING(V,D,RH0P,RH01,RH02,S,XL2,SY1,SY2,THETA, 
+ XN,E1 , CMCRAD.Tl ,T2,WT) 

• 

TIT = TIT + XPV(J) * T1 
T2T = T2T + XPV(J) * T2 
36 CONTINUE 
T1 =T1T 
T2 = T2T 

WT = RH01 * T1 + RH02 * T2 
R12 = CMCRAD 
R22 = CMCRAD + T2 
R1 1 = CMCRAD + T2 + S 
R21 = CMCRAD + T1 +T2 + S 
VB=3.1416*(CMCLEN/1000.)*(R21**2.-R11**2.) 
VW=3.1416*(CMCLEN/1000.)*(R22**2.-R12**2.) 

WTCMC = RHOl * VB + RH02 * VW 

991 CONTINUE 



An Employee-Owned Company 


else IF (NC0DE.EQ.3) then 
C MADDEN 

****** MADDEN MINIMIZES SUM OF THICKNESSES ONLY ******* 


C 45 READ(10,*)D 


c 

c 

c 

c 

c 

c 

c 

c 


C PROJECTILE DIAMETER IN CM **** READ(10,*)D 
READ(22,*,end=l 1)D 
c RE AD( 1 0,*)RHOP 

C ****** (10,*)S 

READ(23,*,end=l 1)S 
c READ(10,*)RHO 
c WRITEO 1,*)’ MADDEN’ 

WRTTE(11 ,*) 

WRITEO 1,*)’ INPUT’ 

WRITE(1 1 *) 

WRITEO L*)’ 

WRITE(11 ,*)’ 

WRITE(11 ,*)’ 

WRTTE(11,*)’ 

WRITE(11 *) 

TIT = 0.0 
T2T = 0.0 


Projectile Diameter In CM = ’ J) 

Projectile Density In IN GM/Cubic CM = ’,RHOP 
Bumper/Wall Density In GM/Cubic CM = ’,RHO 
Bumper/Wall Separation In CM = ’,S 


D0 46 J=l,16 
V = FLOAT(J) 

C CALL MADDEN(V,D,RHOP,S,RHO,Tl ,T2,WT,WTCMC) 

CALL MADDEN(V, D, RHOP, S, RHO, Tl, T2) 

TIT = TIT + T1 * XPV(J) 

T2T = T2T + T2 * XPV(J) 

46 CONTINUE 


c 

c 

c 

c 

c 


Tl =T1T 


T2 = T2T 
WT = Tl +T2 
R12 = 211.0 
R22 = 211.0 +T2 


Rll = 211.0 +T2 + S 


R21 = 211.0 +T1 +T2 + S 
VB=4.27*(R21**2.-R1 1**2.) 
VW=4.27*(R22**2.-R12**2.) 

WTCMC = RHO * (VB + VW) 
c WRITEO 1,*)’ OUTPUT’ 

WRITEO 1,*) 

WRITEO 1,*)’ Bumper Thickness = \T1,’CM’ 
WRITEO 1,*)’ Wall Thickness = \T2,’CM’ 
WRITEO 1,*)’ Minimum Weight = \WT,’CM’ 
WRITEO 1,*)’ CMC Minimum Weight = \WTCM 
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c WRITE(11*) 
c WRITE(11,*) 
c WRITE(11,*) 

C- 

else IF (NC0DE.EQ.4) then 
C WILKINSON 

C **** READ(10,*)D 
READ(22,*,end=l 1)D 
c **** re AD( 1 0,*)RHOP 

READ(22,*,end=l l)RHOP 
C **** (10 ,*)RH01 

READ(24,*,end=l 1)RH01 
C **** (10,*)RHO2 
READ(25,*,end=l 1)RH02 
C **** (10,*)S 
READ(23,*,end=ll)S 
C **** (10,*)XL2 
READ(25,*,end=l 1)XL2 
C **** (10 ,*)CMCLEN 

READ(26,*,end=l 1)CMCLEN 
C **** (10,*)CMCRAD 
READ(26,*,end=l 1)CMCRAD 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


WRITEd!,*)’ WILKINSON’ 


WRITE(11 ,*) 
WRITE(11*)’ 
WRITE(11,*) 
WRITE(11 ,*)’ 
WRITE(11 ,*)’ 
WRITE(11,*)’ 
WRITE(11 ,*)’ 
WRITE(11,*)’ 
WRITE(11,*)’ 
WRITE(11,*) 
TIT = 0.0 
T2T = 0.0 


INPUT’ 

Projectile Diameter In CM = ’ J) 

Projectile Density In GM/Cubic CM = \RHOP 
Bumper Density In GM/Cubic CM = ’.RHOl 
Wall Density In GM/Cubic CM = \RH02 
Bumper/Wall Separation In CM = \S 
Wall Material Constant = ’^CL2 


DO 56 J=l,16 
V = FLOAT(J) 

C CALL WILKINS0N(VJ),RH0P,RH01,RH02,S,XL2, 

C & T1 ,T2,WT,WTCMC) 

CALL WILKINSON(VJ), RHOP, RHOl, RH02.S.XL2, 
& T1.T2) 
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TIT = TIT + T1 * XPV(J) 

T2T = T2T + T2 * XPV(J) 

56 CONTINUE 

T1 =T1T 
T2 = T2T 

WT = RH01 * T1 + RH02 * T2 

R12 = CMCRAD 

R22 = CMCRAD + T2 

Rll = CMCRAD + T2 + S 

R21 = CMCRAD + T1 +T2 + S 

VB=3.1416*(CMCLEN/1000.)*(R21**2.-R11**2.) 

VW=3.1416*(CMCLEN/1000.)*(R22**2.-R12**2.) 

WTCMC = RHOl * VB + RH02 * VW 
c WRITE(11,*)’ OUTPUT’ 
c WRITE(11*) 

c WRTTE(1 1 *)’ Bumper Thickness = ’.Tl.’CM’ 
c WRITE(1 1 *)’ Wall Thickness = ’,T2,’CM’ 
c WRITE( 11*)’ Minimum Weight = ’ ,WT,’GM/Square CM’ 

c WR1TE(1 1,*)’ CMC Minimum Weight = ’.WTCMC, ’KG’ 
c WRITE(11*) 
c WRITE(11*) 
c WRITE(11*) 

C else IF (NCODE.EQ.5) then 
C MODIFIED BURCH 

C **** READ(10,*)D 
READ(22,*,end=l 1)D 
C **** (10,*)RHOl 
READ(24,*,end=l 1)RH01 
C **** (10 ,*)RH02 

. READ(25,*,end=l 1)RH02 
C **** (10,*)S 
READ(23,*,end=l 1)S 
C **** READ(10,*)THETA 
READ(22,*,end= 1 1)THETA 
C **** READ(10 *)XN 
READ(24,*,end=l 1)XN 
C **** (10,*)E1 
READ(24*,end=ll)El 
C **** (10,*)CMCLEN 
READ(26,*,end=l 1)CMCLEN 
C **** (10 ,*)CMCRAD 

READ(26,*,end=l 1)CMCRAD 


***** MODIFIED BURCH ***** 
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c WRITE(11,*)’ MODIFIED BURCH’ 
c WRITEC11,*) 
c WRITE(11,*)’ INPUT’ 
c WRITECll,*) 

c WRITE(11,*)’ Projectile Diameter In CM = ’ J) 
c WRITE(11,*)’ Bumper Density In GM/Cubic CM = \RH01 
c WRITE(1 1,*)’ Bumper/Wall Separation In CM = \S 
c WRITE( 11,*)’ Impact Angle From Normal In Degrees = ’ , THETA 

c WRITE(1 1,*)’ Number Of Plates To Penetrate After First’, ~ ■ 
c + ’ Bumper = ’, XN 

c WRTTE(1 1,*)’ Bumper Youngs Modulus In MSI - ’, El 
c WRITE(11,*) 

El = El * 6.880285E+10 
TIT = 0.0 
T2T = 0.0 

DO 66 J=l,16 

V = FLOAT(J) : 

C CALL BURCH(V,D,RH01,RH02,S, THETA, 

C & XN,E 1 ,T 1 ,T2, WT.WTCMC) 


CALL BURCH(VJD,RH01 ,RH02,S, THETA, 
& XN,E1,T1,T2,T1B,F1) 

TIT = TIT + T1 * XPV(J) 

T2T = T2T + T2 * XPV(J) 

66 CONTINUE 


c 

c 

c 

c 

c 

c 

c 

c 


T1=T1T 
T2 = T2T 

WT = RHOl * T1 + RH02 * T2 

R12 = CMCRAD 

R22 = CMCRAD + T2 

Rll = CMCRAD + T2 + S 

R21 = CMCRAD + T1 + T2 + S 

VB=3.1416*(CMCLEN/1000.)*(R21**2.-R11**2.) 

VW=3.1416*(CMCLEN/1000.)*(R22**2.-R12**2.) 

WTCMC = RHOl * VB + RH02 * VW 
c WRITE(11,*)’ OUTPUT’ 

WRITE(11,*) 

WRITE(1 1 ,*)’ Bumper Thickness = ’,T1,’CM’ 

WRITE(11,*)’ Wall Thickness = ’,T2,’CM’ 

WRITE(1 1,*)’ Minimum Weight = ’.WT.’GM/Square CM’ 
WRITE(1 1 ,*)’ CMC Minimum Weight = ’,WTCMC,’KG’ 

WRITE(11 ,*) 

WRITE(11,*) 

WRITE(11,*) 
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end if 

C HERE WE DEFINE AND CALCULATE NEW VARIABLES 
C NEEDED FOR OUTPUT 

C OPTIMAL BUMPER MASS PER UNIT AREA 
OBMPUA = T1 * RHOl 

C OPTIMAL WALL MASS PER UNIT AREA 
OWMPUA = T2 * RH02 

C OPTIMAL BUMPER RATIO 
OBR = Tl * RHOl /WT 

C OPTIMAL WALL RATIO 
OWR = T2 * RH02 / WT 

C SPACECRAFT INITIAL OPERATING CAPABILITY 
SIOC = I YEAR 1 

C SPACECRAFT MISSION DURATION 
SMD = IYEAR2 - IYEAR1 + 1 


HERE WE WRITE THE CALCULATED OUTPUT VALUES 
TO THE MAIN OUPUT FILE CALLED ’RESULTS.OUT’ 


if ((ibumper_type.eq.l).or.(iwall_type.eq.l)) then 
C Write Calculated Output Values To ’ TABLEl.OUT ’ 

WRITE(37,782) BUMPER_NAME(iiii) ,WALL_NAME(iiii),T 1 , 
+ T2,OBMPUA,OWMPUA,WT,WTCMC 

782 F0RMAT((A14,1X),(A1 1,1X),(5F7.4,1X),F9.2) 


else 


+ 

+ 


+ 

+ 

+ 

+ 


IF (IGRAPH_TYPE.EQ. 1 )THEN 

778 WRITE(27, 779)1, Tl,T2,OBMPUA,OWMPUA,WT,WTCMC,OBR,OWRJI>,RHOP, 
XGROWTHJMONTH 1 JYEAR 1 ,IMONTH2,IYEAR2,ALT,XINCL,XPO,AREAK,S 

WRITE(29,779)I,T1,T2, OBMPUA, OWMPUA,WT,WTCMC,OBR,OWRL>,RHOP, 
XGROWTHJMONTH 1 JYEAR 1 ,IMONTH2,IYEAR2,ALT,XINCL,XPO, ARE AK,S 

779 FORMAT(I9, 1 X, 1 1 (F 1 2.4, 1 X),4(I9, 1X),5(F1 2.4, IX)) 
else if (igraph_type.eq.2) then 

780 WRnE(27,781)I,7,Tl,\\T2,7,OBMPUA,7,OWMPUA,7, 
WT f , f , ,WrCMC, , , , ,OBR, , t , ,OWR,\ , A , , , 3HOP, , > ’, 
XGR0WTH,7JM0NTH1,7JYEAR1,7,IM0NTH2,7JYEAR2,7, 
ALT,’, , ,XINCL, , , , ,XPO, , ,’AREAK,’, , ,S 

WRITE(29,779)I,Tl,T2,OBMPUA,OWMPUA,WT,WTCMC,OBR,OWRJD,RHOP, 
XGROWTHJMONTH 1 .IYEAR 1 .IMONTH2 JYEAR2 ALT,XINCL,XPO,AREAK,S 

781 FORMAT(I9,lX,ll(F12.4,lX),4(I9,lX),5(F12.4,lX)) 
end if 


end if 


1 = 1 + 1 
10 GOTO 20 
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11 CONTINUE 

C HERE WE WRITE THE CALCULATED VALUES OF "V", 
C ”XPV(V)", AND "THETA" TO THE OUTPUT FILE CALLED 
C "PROJECT.OUT" prior to leaving 

WRITE(28,3675)V,XPV(V), THETA 
3675 F0RMAT(3(F9.5,1X)) 


write (*,*) ’ Program Finished’ 

CLOSE (UNIT = 20) 

CLOSE (UNIT = 23) 

CLOSE (UNIT = 26) 

CLOSE (UNIT = 27) 

CLOSE (UNIT = 28) 

CLOSE (UNIT = 29) 

CLOSE (UNIT = 30) 

CLOSE (UNTT = 35) 

CLOSE (UNIT = 37) 

STOP 

END 

C 

C C 

C SUBROUTINES BEGIN HERE 

C 

SUBROUTINE DEBRIS(XGROWTH,SOLAR,XPSI,IMONTHl,IYEARl,IMONTH2, 
& IYEAR2ALT,XINCL^P0,AREAKD,XPV,IVMAX) 

C DIMENSION SOLAR(100),XPV(100)^CPSIV(105) <-- MODIFIED 
DIMENSION SOLAR(l 188), XPV(IOO) 

G1TOT = 0.0 
G2TOT = 0.0 

NYEARS = IYEAR2 - 1 YEAR 1 + 1 

NMONTHS = 12* (IYE AR2-IYEAR 1 )+imonth2-imonth 1 

DO 582 EL=1, NMONTHS 
XPHI 1 = 10.* *(( ALT/200.)- 

& (SOLAR( 1 2*(IYEAR 1-1987-1 )+imonth 1 +IL- 1 )/140.)- 1 .5) 

XPHI = XPHI1 / (XPHI1 + 1.0 ) 

Gl=(l.+2.*XGROWTH)**(IYEARl-1985+(imonthl+il-2)/12.0) 

G2=(l.+XGROWTH)**(TYEARl-1985+(imonthl+il-2)/12.0) 

G1TOT = G1TOT + XPHI * G1 
G2TOT = G2TOT + XPHI * G2 
582 CONTINUE 
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FLUX = 12.0 * ALOG(XPO) / (AREAK * XPSI) 

DEN = -1.0 * (5.9499E-07 * G2TOT + FLUX) 

XNUM = .0000105 * G1TOT 
D=(XNUM/DEN)**0.4 
YG - 250.0 
YF = 0.0 
YC = .0125 

YE = .55 + .005 * (XINCL - 30.0 ) 

YH=1.0-0.0000757*(XINCL-60.0)**2.0 

YA - 2.5 

YB-.3 

YD = 1.3 - .01 * (XINCL - 30.0 ) 

YVO = 7.7 

EF(XINCL .LE. 60.0)THEN 
YB = .5 
YG = 18.7 

YVO = 7.25 + .015 * (XINCL - 30.0 ) 

END IF 

IF(XINCL .LE. 80.0 .AND. XINCL .GT. 60.0)THEN 
YB = .5 - .01 * (XINCL - 60.0 ) 

YG=18.7+0.0289*(XINCL-60.0)**3.0 
END IF 

IF(XINCL .GT. 100.0)THEN 
YC = .0125 + .00125 * (XINCL - 100.0 ) 

END IF 

IF(XINCL .LE. 50.0)THEN 
YF=0.3+0.0008*(XINCL-50.0)**2.0 
END IF 

IF(XINCL .GT. 50.0 AND. XINCL .LE. 80.0)THEN 
YF = .3 - .01 * (XINCL - 50.0 ) 

END IF 

XSUMIV = 0.0 
IVMAX = 1 
IV = 1 

584 XPV(TV)=YG*2.7183**(-1.0*((IV-YA*YV0)/(YB*YV0))**2.0) 

XPV(IV)=XPV(rV)+YF*2.7183**(-1.0*((rV-YD*YV0)/(YE*YV0))**2.0) 
XPV(W)=XPV(IV)*(2.0*IV*YV0-IV**2.0) 
XPV(IV)=XPV(IV)+YH*YC*(4.0*IV*YV0-IV**2.0) 

IF(XPV(TV) .LE. 0.000)THEN 
XPV(IV) = 0.0 
IVMAX = IV 
GOTO 586 
END IF 

XSUMIV = XSUMIV + XPV(IV) 

IV = IV + 1 
GOTO 584 

586 DO 588 1=1, IVMAX 
XPV(I) = XPV(I) / XSUMIV 
588 CONTINUE 
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RETURN 


END 

C 

SUBROUTINE DEBRISOLD(T, XPO, AREAK, D) 

FLUX = -1.0 * ALOG(XPO) / (AREAK * T) 

F = ALOG 10(FLUX) 

C***** MS-FORTRAN DOES NOT ALLOW CONSECUTIVE MATHEMATICAL 
C OPERATORS TO BE PLACED ADJACENT TO ONE ANOTHER, i.e. 

C YOU CAN *** NOT *** HAVE: 

C 

c 

C IF(F.GE.-5.46)THEN 

C D= 1 0. * * ((F+5 .46)/-2.52) 

C END IF 

C IF(F.GE.-5.9.AND.F.LT.-5.46)THEN 

C D=10.**((F+5.02)/-0.44) 

C END IF 

C IF(F.LT.-5.9.ANDE.GE.-6.4)THEN 

C D=10.**((F+5.78)/-0.063) 

C END IF 

C IF(F.GE.-7.0.AND.F.LT.-6.4)THEN 

C D=10.**((F+6.33)/-0.0067) 

C END IF 

C IF(F.GE.-7.3.AND.F.LT.-7.0)THEN 

C D=10.**((F+6.88)/-0.0012) 

C END IF 

C IF(F.GE.-7.6.AND.F.LT.-7.3)THEN 

C D=10.**((F+6.6)/-0.002) 

C END IF 

C IF(F.GE.-8.0.AND.F.LT.-7.6)THEN 

C D=10.**((F+5.6)/-0.004) 

C END IF 

C 

IF (F.GE.- 5.46) THEN 
D=10.**((F+5.46)/(-2.52)) 

END IF ^ ' ‘ = 

IF(F.GE.-579AND.FI.T.-5.46)THEN 

D=10.**((F+5.02)/(-0.44)) 

END IF 

IF(F.LT.-5.9.AND.F.GE.-6.4)THEN 

D=10.**((F+5.78)/(-0.063)) 

END IF 

IF(F.GE.-7.0.AND.F.LT.-6.4)THEN 

D=10.**((F+6.33)/(-0.0067)) 

END IF 

IF(F.GE.-7.3.AND.F.LT.-7.0)THEN 

D=10.**((F+6.88)/(-0.0012)) 

END IF 
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IF(F.GE.-7.6.AND.F.LT.-7.3)THEN 

D=10.**((F+6.6)/(-0.002)) 

END IF 

IF(F.GE.-8.0.AND.F.LT.-7.6)THEN 

D=10.**((F+5.6)/(-0.004)) 

END IF 

RETURN 

END 

C 

C SUBROUTINE METEOROID(SA,T,PO,ALT,DENS,D) 

SUBROUTINE METEOROID(AREAK, T, XPO, ALT, DENS, D, L) 


T = 31536000.0 * T 

FLUX = -1.0 * ALOG(XPO) / (AREAK * T) 

RA = 6371.0 / (6371.0 + ALT) 

GE * .568 + .432 * RA 

THETA = ATAN(6371.0 / SQRT(ALT * (ALT + 2.0 * 6371.0 ))) 

S = (1-0 + COS (THETA)) / 2.0 
FLUX = FLUX / (GE * S) 

F = ALOG10(FLUX) 

IF (F.GE. - 4.403) THEN 

c WRITE(1 1,*)’ MASS IS TOO SMALL’ 

GOTO 1001 
END IF 

IF(F.GT.-7. 103.AND.F.LT.-4.403)THEN 
RAD = 2.509 - .25 * (14.339 + L) 

XM=10.**((-1.584+SQRT(RAD))/.125) 

END IF 

IF(F.LE.-7. 103. AND.F.GE.- 14.37)THEN 
C XM=10.**((14.37+F)/-1.213) <— NOT ALLOWED IN MS-FORTRAN 
XM=10.**( (14.37+F)/(-1.213) ) 

END IF 

IF (F.LT. - 14.37) THEN 

c WRITE(1 1 *)’ MASS IS TOO LARGE’ 

GOTO 1001 
END IF 

D=(1.91*XM/DENS)**.333 
1001 CONTINUE 
RETURN 

END 

C SUBROUTINE NYSMITH(V,D,H,RH01,RH02,T1,T2,WT,WTCMC) 
SUBROUTINE NYSMITH(V, D, H, RHOl, RH02, Tl, T2) 

C DMAX=0.24*H*V**-0.2 <— NOT ALLOWED IN MS-FORTRAN 
DMAX=0.24*H*V**(-0.2) 
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IF (D.GT.DMAX) THEN 

WRTTEQ 1 ,*)’ NO SOLUTION-PROJ. DIA. TOO LARGE FOR NYSMITH’ 

else 

T1-(1.93*V**0.18*D**1.91/H**0.91)*((RH02/RH01)**0.65) 

T2 = 1.86 * T1 * RHOl / RH02 
END IF 

RETURN 

END 


***** PEN4 ***** 

C SUBROUTINE BOEING(V,D,RHOP,RHO 1 ,RH02,S,XL2,SY 1 ,SY2, THETA, 

C & XN,E1,CMCRAD,T1,T2,WT,WTCMC) 

SUBROUTINE BOEING(V ,D,RH0P,RHOl,RHO2,S,XL2,S Y1 ,S Y2,THETA, 

+ XN,E1,CMCRAD,T1,T2,WT) 

T1 = .16 
V-V* 3280.0 
D = D/ 30.48 - • 

RP-D/2.0 
RHOP = RHOP* 1.94 
RHOl = RHOl * 1.94 
RH02 = RH02 * 1.94 
NITSP = 0 
NITSP = NITSP + 1 
NP1 = 0 

TIP = T1/ 30.48 

T2P * FT2P(RHOP,V, RP, SY1, THETA, RH01,SY2,D,RH02,T1P) 

WT = RHOl * TIP + RH02 * T2P 
IF (NITSP.EQ.1) THEN 
TIPI = 1.1 * TIP 

T2P1 = FT2P(RHOP,V,RP,SYl, THETA, RHOl, SY2JD,RH02,T1P1) 

WT1 = RHOl * TIPI + RH02 * T2P1 
END IF 

IF (WT1.GT.WT) THEN 
TIPI = .82 * TIPI — - 

T2P1 = FT2P(RH0P,V,RP,SY1,THETA,RH01,SY2J),RH02,T1P1) 

WT1 = RHOl * TIPI + RH02 * T2P1 
590 IF (WT1 .GT.WT) THEN 
GOTO 601 

ELSE 

TIP = TIPI 
T2P = T2P1 
WT = WT1 
NP1 = NP1 + 1 

IF (NP1.EQ.100) THEN 

c*** WRITE(11 ,*)’ NO CONVERGENCE IN PEN4’ ' 

GOTO 557 
END IF 
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TIPI = .9 * TIPI 

T2P1 = FT2P(RHOP,V,RP,S Y 1 ,THETA,RHO 1,SY2 J),RH02,T1P1) 
WT1 = RH01 * TIPI + RH02 * T2P1 
GOTO 590 
END IF 
FI <?F 

579 TIP = TIPI 
T2P = T2P1 
WT = WT1 
TIPI * 1.1 * TIPI 

T2P1 = FT2P(RHOP,V,RP,SYl, THETA, RH01,SY2D,RH02,T1P1) 
WT1 = RHOl * TIPI + RH02 * T2P1 
IF (WT1.GT.WT) THEN 
GOTO 601 
FT <?F 

NP1 = NP1 + 1 

IF (NP1.EQ.100) THEN 

c*** WRITE(11*)’ NO CONVERGENCE IN PEN4’ 

GOTO 557 
END IF 

GOTO 579 
END IF 
END IF 

601 CONTINUE 
D = 30.48 * D 
RHOP = RHOP / 1 .94 
RHOl = RHOl / 1.94 
RH02 = RH02 / 1 .94 
TIP = 30.48 * TIP 
T2P = 30.48 *T2P 
IF(T1P/D.LE.0.4)VF=4100 
IF(T1P/D.GT.0.4)VF=4986*(T1P/D)**0.21 
VF = VF + 4000.0 

IF (V.LE.VF) THEN 

c*** WRITE(1 1 *)’ INSIDE OF PEN4 LIMITS’ 

T1 =T1P 
T2 = T2P 
GOTO 1102 
END IF 

557 CONTINUE 

***** WILKINSON ***** 

V = V/ 3280.0 

Tl=0.604*D**2.*RHOP/(S*RHOl) 

T1 =T1 * SQRT(V * COS(THETA) / XL2) 

T2 = T1 * RHOl / RH02 
RATIO = D * RHOP/ (T1 * RHOl) 

IF (RATIO.GT.1.0) GOTO 1458 
IF(RATIO.LE.1.0)T2=T2/RATIO 
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1458 CONTINUE 

c*** WRrrE(ll,*)’TlW = \T1,’T2W = \T2 
***** MODIFIED BURCH ***** 

THI=0.8 16*(0.5236*RHOP*D**3.0)**0.352*(RHOP**0. 167) 
THI=THI*(V**0.875)/(0.8467*RH01**0.5) 

THI = THI/2.54 
TLO = 0.0 

TINTERVAL = THI - TLO 
C 

ISEED=9 14 1 1 <-- THIS IS NO LONG ER N EEDED, SINCE 
WE WILL USE OUR OWN SEED VALUE 
WHICH WILL COMPLY WITH MS-FORTRAN 

VB = V * 3280.0 
DB = D/2.54 
CM = SQRT(E1 /RHOl) 

CM = CM/ 30.48 
SB = S / 2.54 
RHOP = RHOP * .036215 
RHOl = RHOl * .036215 
RH02 = RH02 * .036215 

IF (THETA.LE.0.001) GOTO 125 
CHI = TAN(THETA) - .5 
XPENALTY =1.0 
TIB = THI 

F1=2.42*(DB/T1B)**0.333+4.26*(T1B/DB)**0.333-4.18 
F2 = .5 - 1.87 * (TIB / DB) + (5.0 *T1B/DB- 1.6) 

+ * CHI * CHI * CHI 
F2 = F2 + (1.7 - 12.0 * TIB / DB) * CHI 
F3=0.32*(T1B/DB)**0.83 

F3=F3+0.48*(T1B/DB)**0.33*(SIN(THETA))**3.0 
c*** WRITE(11,*)’DB = ’£>B,’XN = ’^CN.’Fl = ’,Fi,’F2 = 
c*** WRITE(1 1,*)’THI = THI, ’CHI = ’.CHI 
IF (FI + .63 * F2.LT.0.001) THEN 
T2F = 2116.8 * CMCRAD / SY2 
GOTO 483 
END IF 

T2F=D B*((F1+ 0.63*F2)/XN)** 1 .7143 
T2F=T2F*(CM/VB)**2.2857 
T2F=T2F*(DB/SB)**0.7 143 
483 XNN=F3*(DB/T2F)*(CM/VB)** 1 .333 

YDELTA = 0.0 

IF(XNN.GT.0.850)YDELTA=1 .000 
TOTPEN=YDELTA*XPENALTY*(XNN-0.85)**2.00 

WTB = RHOl * TIB + RH02 * T2F + TOTPEN 

WTMIN = WTB 
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T1BEST = THI 
T2BEST = T2F 
TOTPENBEST = TOTTEN 
482 DO 48 1 IPENALTY= 1 ,460 
C 

C T1 B=TINTER VAL*RAN(ISEED) <-- RAN IS NOT USED IN 

CMS -FORTRAN 

CALL RANDOM(RANVAL) 

TIB = TINTERVAL * RANVAL 
F 1 =2.42* (DB/T 1 B) * *0.333+4.26* (T 1 B/DB)**0.333-4. 1 8 
F2 = .5 - 1.87 * (TIB / DB) + (5.0 * TIB / DB - 1.6) 

+ * CHI* CHI* CHI 
F2 = F2 + (1.7 - 12.0 * TIB /DB) * CHI 
F3=0.32*(T 1 B/DB) * *0. 83 

F3=F3+0.48*(T1B/DB)**0.33*(SIN(THETA))**3.0 
IF (FI + .63 * F2.LT.0.001) THEN 
T2F - 21 16.8 * CMCRAD / SY2 
GOTO 484 
END IF 

T2F=DB*((F1+0.63*F2)/XN)** 1.7143 
T2F=T2F*(CM/VB)**2.2857 
T2F=T2F*(DB/SB)**0.7143 
484 XNN=F3*(DB/T2F)*(CM/VB)** 1.333 

YDELTA * 0.0 

IF(XNN.GT.0.850)YDELTA=1.000 
TOTPEN = YDELTA*XPENALTY* (XNN -0. 85)**2.00 
WTB = RHOl * TIB + RH02 * T2F + TOTPEN 
IF (WTB.LT:WTMIN) THEN 
WTMIN = WTB 
T1BEST = TIB 
T2BEST = T2F 
TOTPENBEST = TOTPEN 
END IF 

481 CONTINUE 

IF (TOTPENBEST.GT.O.OO 1 ) THEN 
XPENALTY = XPENALTY * 10.0 
IF (XPENALTY.GT.1.0E12) THEN 
GOTO 485 
END IF 
GOTO 482 
END IF 

485 TIB = T1BEST 
T2B = T2BEST 

c*** WRITE(1 1 ,*)’T1B = \T1B,’T2B = ’,T2B,’K = ’^PENALTY 

c ***+* WRITE(1 1*) ’TOTPENBEST = TOTPENBEST 

GOTO 499 
125 CONTINUE 
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WRITE(11 *)’RH01 = \RH01 
c***** WRITE(1 1 ,*)’RH02 = ’,RH02 


XK1=(DB/XN)** 1 .7 1 *(CM/VB)**2.29/SB**0.7 1 
VDELTA = 0.0 
DELTA3 = .52 

1099 DELTA2 = 2.33 * (1.0 - 1.57 * DELTA3) 

DELTA 1 = 1.33 * (2.0 * DELTA3 - 1.0 ) 

VDELTA1=(1VDELTA1)**DELTA1*(2.8*XK1/(DELTA2*DB**0.57)) 

+ **DELTA2 

VDELTA 1 =VDELTA1 *( 1 .58*XK1 *DB**0.57/DELTA3)**DELTA3 
VDELTA 1 =VDELTA 1 * (RHO 1 **DELTA 1 )*(RH02**(DELTA2+DELTA3)) 
IF (VDELTA1.LT. VDELTA) THEN 
DELTA 1 = 1.33 * (2.0 * DELTA3 - 1.04) 

TIB = DELTA 1 * VDELTA / RHO 1 — • — 

T2B = (VDELTA - TIB * RHOl) / RH02 
GOTO 499 
END IF 

VDELTA = VDELTA 1 
DELTA3 = DELTA3 + .02 

IF (DELTA3.GT.0.63) THEN 
TIB = DELTA 1 * VDELTA / RHOl 
T2B = (VDELTA - TIB * RHOl) / RH02 
GOTO 499 
END IF 

GOTO 1099 
499 CONTINUE 

***** COMPARISON OF MODIFIED BURCH AND WILKINSON ***** 
199 CONTINUE 
T10W = T1 /2.54 
IF (THETA.LT.0.00 1 ) GOTO 486 

F10W=2.42*(DB/T10W)**0.333+4.26*(T10W/DB)**0.333-4.18 
F20W = .5 - 1.87 * (T10W/DB) + (5.0 * T10W/DB - 1.6) 

+ * CHI * CHI * CHI 

F20W = F20W + (1.7 - 12.0 * T10W / DB) * CHI 
F30W=0.32*(T10W/DB)**0.83 

F30W=F30W+0.48*(T10W/DB)**0.33*(SIN(THETA))**3.0 
IF (F10W + .63 * F20W.LT.0.001) THEN 
T2FT10W = 21 16.8 * CMCRAD / SY2 
GOTO 487 
END IF 

T2FT10W=DB*((F10W-K).63*F20W)/XN)**1.7143 
T2FT10W=T2FT10W*(CMyVB)**2.2857 
T2FT10W=T2FT10W*(DB/SB)**0.7 143 
487 T2BT10W = T2FT10W * 2.54 

XNNT 1 0W=F30W* (DB/T2FT 1 0W) * (CM/VB)* * 1 .333 
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IF (XNNT10W.GT.0.85) THEN 
T2BT10W = 0.0 
END IF 

RATIOB = (DB * RHOP) / (RHOl * TIB) 
T2WT10B=0.364*D**3.*RHOP*V*COS(THETA)/(XL2*RHO2*S**2.) 
IF(RATIOB.GT. 1.0)T2WT10B=T2WT10B*RATIOB 
IF(T2BT 1 0W.GT.T2)T2=T2BT 1 0W 
T2B = T2B * 2.54 

EF(T2WT 1 0B.GT.T2B)T2B =T2WT 1 OB 
TIB = TIB *2.54 
RHOP = RHOP/. 036215 
RHOl = RHOl / .036215 
RH02 = RH02 / .036215 

IF (RHOl * TIB + RH02 * T2B.LT.RH01 * T1 + RH02 * T2) THEN 
T1 = TIB 
T2 = T2B 
END IF 

GOTO 155 

486 F10W=1.58*(DB/T10W)**0.57+2.80*(T10W/DB)**0.57 

T2BT10W=(F10W/XN)** 1 .7 1 *(CM/VB)**2.29*DB** 1 .71 
T2BT10W=T2BT10W/SB**0.7 1 
T2BT10W = T2BT10W * 2.54 
RATIOB = (DB * RHOP) / (RHOl * TIB) 

T2WT10B=0.364*D**3.*RHOP*V*COS(THETA)/(XL2*RHO2*S**2.) 
IF(RATIOB.GT. 1 .0)T2WT10B=T2WT10B*RATIOB 
EF(T2BT 1 0W.GT.T2)T2=T2BT 1 0W 
T2B = T2B * 2.54 

IF(T2WT10B.GT.T2B)T2B=T2WT10B 
TIB = TIB * 2.54 
RHOP = RHOP/. 036215 
RHOl = RHOl / .036215 
RH02 = RH02/. 036215 

IF (RHOl * TIB + RH02 * T2B.LT.RH01 * T1 + RH02 * T2) THEN 
T1=T1B 
T2 = T2B 
END IF 

155 CONTINUE 

1 102 IF (T2.LE.0.01) THEN 

T2 = 21 16.8 * CMCRAD / SY2 
c***** WRITE(1 1,*)’T1P = \T1,’T2P = ’,T2 
END IF 

c*** WRITE(1 1 *)’T1 = ’,T1,’T2 = \T2 

156 RETURN 
END 

C 

FUNCTION FT2B (DB, TIB, XN, CM, VB, SB) 
F1=2.42*(DB/T1B)**0.33+4.26*(T1B/DB)**0.33 
FI = FI -4.18 

FT2B=(Fl/XW**1.71*(CM/VB^**2.29*DB**1.71/SB**njhi m Ml" 
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RETURN 


END 

C 

FUNCTION FT2P (RH0P,V,RP ) SY1,THETA,RH01,SY2,D,RH02,T1P) 

A=1.33*RHOP*(V*RP)**2. ..... . 

B = 8.0 * SY1 * EXP(- .0003 125 * V) / COS(THETA) 

C=1.33*RHOP*RP**2.0 

D1 = RP * RHOl / COS (THETA) 

XK1=1.67*(RHOP/(2.*SY2))**0.31 
XKl=XKl*(0.281*D*RHOP/RHO2)**0.33 
XK1 = XK1 * COS(THETA) 

C1P1 = (A - B * TIP) / (C + D1 * TIP) 

IF (C1P1.LE.0.001) THEN 
FT2P = 0.0 
GOTO 999 
END IF 

FT2P=XK1*C1P1**0.31 
999 RETURN 

end . .. .. .. .. ,,,,, .. ... . . ' 

c 

C SUBROUTINE MADDEN(V,D,RHOP,S,RHO,Tl,T2,WT,WTCMC) 
SUBROUTINE MADDEN(V, D, RHOP, S, RHO, Tl, T2) 


V = V * 100000.0 

T1=0.009*SQRT(V)*RHOP*D**2.0 

Tl=Tl/(S*RHO**1.5) 

T2 = Tl 

RETURN 

END 

C 

C SUBROUTINE WILKINS0N(V > D ) RH0P ) RH01,RH02,S > XL2 ) 
C & Tl ,T2,WT,WTCMC) 

SUBROUTINE WILKINS0N(VX>,RH0P,RH01,RH02,S^L2, 
& T1.T2) 


T1^.604*D**2.*RHOP/(S*RHOl) 
Tl = Tl * SQRT(V / XL2) 

T2 = Tl * RHOl / RH02 
RATIO = D * RHOP/ (Tl * RHOl) 

IF (RATIO.GT.1.0) GOTO 3683 
IF(RATIO.LE. 1 .0)T2=T2/RATIO 
3683 CONTINUE 
RETURN 

END 


C 

***** MODIFIED BURCH ***** 
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C SUBROUTINE BURCH(VJ),RH01,RH02,S,THETA, 

C & XN,E 1 ,T1 ,T2, WT, WT CMC) 

SUBROUTINE BURCH(VJD,RH01,RH02,S,THETA, 

& XN,E1,T1,T2,T1B,F1) 

VB = V * 3280.0 
DB = D / 2.54 
CM = SQRT(E1 / RHOl) 

CM = CM/ 30.48 
SB = S/2.54 

IF (THETA.LE.0.00 1 ) GOTO 425 
CHI = TAN(THETA) - .5 
F2=0.5-1.87*(T1B/D)+(5.*T1B/D-1.6)*CHI**3.0 
F2 = F2 + (1.7- 12.0 * TIB /D) * CHI 
F3=0.32*(T1B/D)**0.83 

F3=F3+0.48*(T1B/D)**0.33*(SIN(THETA))**3.0 
T2F=D*((F1+0.63*F2)/XN)*(CM/V)**2.29 
T2F=T2F*(D/S)**0.71 
T2N=F3*(CM/V)** 1 .33*D/XN 
IF(T2N.GE.T2F)T2B=T2N 
IF(T2N.LT.T2F)T2B=T2F 
T2B = T2B * 2.54 
IF(T2B.GT.T2)NREGION=3 
IF(T2B.GT.T2)T2=T2B 
GOTO 499 
425 CONTINUE 
NITSB = 0 

XK1 =(DB/XN)** 1 .7 1 *(CM/VB)**2.29/SB**0.7 1 
VDELTA = 0.0 
DELTA3 = .52 

1099 DELTA2=2.33*(1.-1.57*DELTA3) 

DELTA 1 = 1.33 * (2.0 * DELTA3 - 1.0 ) 

VDELTA 1 =( 1 ./DELTA 1 )* *DELTA 1 * (2. 8 *XK 1/(DELTA2*DB * *0.57)) 

+ **DELTA2 

VDELTA 1 =VDELTA 1 *(1.58*XK1 *DB**0.57/DELTA3)**DELTA3 
VDELTA 1 =VDELTA 1 *(RHO 1 * *DELTA 1 ) *(RH02**(DELTA2+DELTA3)) 
IF (VDELTA 1 .LT. VDELTA) THEN 
DELTA 1 = 1.33 * (2.0 * DELTA3 - 1.04) 

T1 = DELTA 1 * VDELTA / RHOl 
T2 = (VDELTA - T1 * RHOl) / RH02 
GOTO 499 
END IF 

VDELTA = VDELTA 1 
DELTA3 = DELTA3 + .02 
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184 


v. 


IF (DELTA3.GT.0.63) THEN 
T1 = DELTA 1 * VDELTA / RH01 
T2 * (VDELTA - T1 * RHOl) / RH02 
GOTO 499 
END IF 

GOTO 1099 
499 CONTINUE 
T1=T1 *2.54 
T2 = T2 * 2.54 
RETURN 
END 
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